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ABSTRACT 

Inversion of Laplace transforms has been accomplished 
by a numerical integration along appropriate paths in the 
complex plane. Two general procedures have been used. 

The simpler and more economical employs a simple path, such 
as a parabola, which bends to the left. Accuracy is main- 
tained by monitoring the oscillation of the integrand. A 
second method employs a steepest descent contour. 
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INTRODUCTION 



I . 

If F(t) is an integrable function of the real variable 
t, the function f(s) defined by the integral 

e -St F(t)dt (1) 

is called the Laplace transformation, or the Laplace trans- 
form, of F(t) and is frequently indicated by the notation 

f (s) =£[F(t)] (2) 

The Laplace transformation is a linear integral trans- 
formation which is widely used in applied mathematics and 
technology. A typical application is as follows. An un- 
known function F(t) satisfies a certain differential equa- 
tion and specified end conditions. Employing theorems 
applicable to the Laplace transform, a function f(s) is 
determined as the solution of an algebraic equation which 
corresponds, in an appropriate fashion, to the differential 
equation and end conditions of F(t). Then the problem is 
reduced to the following task: having established f(s), 

it is necessary to determine the function F(t) for which 
f (s) -&[F(t)] . 
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This problem is called finding the inverse (Laplace) 



transform of the function f(s). This is customarily 



written as 



F(t) =<£ 1 [f(s)] 



(3) 



It may be shown that the following inversion integral 



accomplishes the desired inversion. -In this expression 
s is the complex variable 



and the integration is along a suitable path in the complex 
s-plane . 

If in fact there is a function F(t) such that the given 
f (s) = &[F(t)] , then the function determined by equation 
(4) is, indeed, the function F(t). The result is unique 
in the following sense. If 



F(t) 




f ( s ) ds 



(4) 



c 



s = x + iy 



(5) 



CO 



f (S) 




e 



-st 




( 6 ) 



o 



and 




c 
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then F^t) and FgCt) differ only on a set of Lebesgue 
measure zero. 

In this thesis numerical methods will be utilized to 
perform the integration of equation (4) and functions for 
which the inversion is successful will correspond to 
functions F(t) which are at least piecewise continuous, 
having only isolated discontinuities, if any are present 
at all. If F(t) is discontinuous for t = a, then the 
recovered function satisfies the condition 

F( a) = Lim[F( a+e ) + F(a-e)] (8) 

e -*■ 0 

The path of integration, indicated by C in equation 
(4) is frequently described as an infinitely long, vertical 
straight line x = constant, to the right of any singulari- 
ties of f(s). (A theorem of complex variables shows that 
f(s) is analytic in its region of convergence, a right-hand 
plane, except for isolated singularities.) This straight 
vertical contour is frequently called the Bromwich contour. 

Usually inversion of a Laplace transform, that is find- 
ing F(t) =«£~l[f(s)], is accomplished by use of a table of 
transform pairs. One of the most complete of these is that 
of Roberts and Kaufman [14] . Other such tables are listed 
as references [5], [6], [7], and [13]. 

Several theorems concerning the Laplace transform extend 
the usefulness of such transformation pair tables. For 
example 
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1 [f 1 (s) + f 2 ( s ) ] 1 ff 1 (s)] 1 [f 2 ( s )] 

( 9 ) 

However, it is not an infrequent occurrence that one is 
called upon to find F(t) in a case where the tables and 
the theorems are of no direct help. There has developed 
a considerable literature concerning this problem. 

One obvious method is to attempt to expand the given 
f(s) as a series of functions f^(s)for which the inverses 
are known. Then using the linear property indicated in 
equation (9), the result may be expressed as 

c£ _1 [f(s)] =d£" 1 [Ea i f i (s)] = ZaX _1 [f.(s)] (10) 

The integration indicated in equation (4) can also be 
performed by various analytical and approximate methods. 

In particular, one can perform the integration along the 
Bromwich contour by a numerical procedure. Much of the 
literature is devoted to such methods. 

In general these methods amount to the following. A 
finite number of points along the Bromwich contour are 
determined according to some law and the integrand is 
evaluated at such points. There are associated weight fac- 
tors and the integral is approximated as a sum of products 
of integrand values times weight factors. One of the most 
widely used of such methods is described by Salzer [15], 
[16], and [17]. Other procedures of the type are treated 
in [18], [20], [2], [8], [9], [10], and [21]. 
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It is interesting to note that the same classical poly- 
nomials which are useful in the expansions indicated by 
equation (10) are also encountered in locating points for 
evaluation using a numerical integration scheme so that 
the final algorithm may be the same even though the original 
motivation was quite different. 

Other methods which do not fall clearly into the 
classes of methods discussed above have been presented by 
Widder [22], and Bellman, Kalaba, and Lockett [3]. 

These approximate methods have not always been success- 
ful. Hiep [12] employed Salzer's method for numerical in- 
version on a problem of heat transfer in porous media and 
found it difficult to obtain results which did not exhibit 
physically incorrect behavior, e.g., in the vicinity of a 
point of engineering interest, the temperature of the 
medium fell below sink temperature and rose above source 
temperature. In similar work on a conjugated heat transfer 
problem Zargary [23] encountered the same difficulty, and, 
after devoting considerable time and attention to the prob- 
lem, was forced to abandon the use of numerical inversion. 

To overcome such difficulties in problems of these 
kinds, we have reverted to direct numerical integration of 
equation (4) along appropriate contours in the s-plane. 

Our procedures, which appear thus far to be efficient and 
reliable, are described in what follows. 
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II. DISTORTED CONTOURS 



We have found that great advantage accrues from distort- 
ing the Bromwich contour in ways which will subsequently 
be described. 

According to theorems concerning analytic functions of 
a complex variable, the Bromwich contour can be distorted 
arbitrarily as long as the process of distortion does not 
cause the contour C to cross over a singularity of f(s) and 
as long as the two ends of C extend to infinity in a manner 
which preserves the convergence of the integral in equation 
( 4 ). 

Thus, with the supposition that f(s) is analytic, with 
the exception of isolated singularities, and is real valued 
for real s, let g(s) represent the integrand 



g ( s ) = e St f ( s ) = u(s) + iv( s ) 



(ID 



where u(s) and v(s) are the real and imaginary harmonic 

components of g(s). Furthermore, assume that the contour C 

* * * 



consists of 
in Figure 1. 
over arcs ds 



the two symmetrical parts C and C as shown 
Consideration of the contributions to £ 1 [ f (s) ] 



and ds 9 


on C 


dx 


= -dx 


1 CM 


= dy 


u 0 


= u. 


2 


1 


v 2 


= " v l 



* * 

1 

1 



( 12 ) 
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FIGURE 1. Distorted Bromwich Contour. Contour is 

distorted into C* and its reflection C**. 



Thus , 



1 

2iri 



[gl ds l 



+ g 2 ds 2 ] “ 



dx. 



+ u 1 dy 1 ] 



(13) 



and, therefore, 



F(t) = ^ J [ vdx + udy ] (14) 

* 

C 

* 

Equation (14) is, of course, further simplified if C 
is selected to be a path along which v=0. Additional 
theoretical reasons exist for such a selection. Such a path 
will subsequently be referred to as a path of steepest 
descent. Numerical integration along such a path, or a 
close facsimile thereof, was done by Esch [11] in 1957 and 
Carrier, Krook, and Pearson [4] recommended further 
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exploration of this procedure. However, no further such in- 
vestigation seems to have been done until the present thesis. 

In Chapter V, a description of an algorithm for performing 

* 

numerical integration along a path C , for which v=0 over 

a substantial portion of its range, will be presented. 

Almost by accidental discovery, the present work has 

revealed that it is computationally more efficient simply 

* 

to make an arbitrary choice for the contour C which satis- 
fies the following criteria: 

* 

(1) C is a smooth continuous curve of the form 

s ( p ) = X(p) + iy(p) (15) 

where p is a real parameter 0 <_ p _< °°. 

(2) C lies to the right of all singularities of f(s). 

(3) C does not extend sufficiently far to the right 

to cause computational difficulties due to the 

factor e St . 

* 

(4) C approaches infinity in such a way that x 
approaches -°° and y approaches 00 . 

(5) Oscillation of the integrand is not excessive. 
(Otherwise there may be loss of accuracy due to 
the sampling rate of the integration algorithm 
seriously mismatching variations in the integrand 
and/or positive and negative contributions nearly 
cancelling each other.) 

Such a contour will be referred to as a simple parameterized 
contour. Discussion of how to assure satisfying criteria 
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(2), (3), and (5) will be deferred until after the develop- 
ment of the actual algorithm. 

The integral shown in equation (14) is approximated by 
the simple trapezoidal sum 

M 

F(t) = 23 {[v(s k )+v(s k + l ))[x k + l - x kl 

k=l 

+ [u(s k )+u(s k+1 )][y k+1 - y k ]} (16) 

where 

s k = x k + iy k = x(p k )+ ly(p k> (17) 

The arbitrarily chosen functions x(p) and y(p) are monotoni 

cally increasing and p. , p„,..., p is an increasing 

1 z n 

sequence . 

A very simple example, but one which has been employed 
very successfully, is given by the following equations 

x = A - Bp 2 

y = P (18) 

p = 0, A, 2A , . . . 

The real number A is selected so as to satisfy criteria (2) 
and (3) above. The real number B has, usually, been 
assigned a value of one. 

The sum includes only a finite number of terms. It is 
necessary to provide an appropriate criterion for terminat- 
ing the process of summation. This will be discussed later 
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In order to monitor the possibility of inaccurate 
evaluations resulting from oscillations, in addition to the 
sum F(t) given by equation (16), the sum G(t) given by 



G(t) = 



_ 1 _ 

2tt 



M 



k=l 



{(v(s k+ i)-v(s k )][x k+ i-x k ] 



+ t u(s k + l )+u(s k )]ly k + l- y k 1) 



(19) 



is also obtained. This is the sum of the absolute values 
of the addends to the sum (16). Additionally, the number 
of times successive addends are of opposite signs is 
recorded. 

The present numerical experiments, to be described 

later, were generally so successful that the effect of 

oscillation could not be seen. However, in some cases 

* 

where poor contours C were deliberately chosen, it was 
found that inaccurate results were obtained if G(t) was 

g 

many orders of magnitude (e.g., 10 ) times as great as F(t). 
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Ill . 



T HE COMPUTER PROGRAM AND 
ITS PRINCIPAL SUBROUTINES 



Previous discussion proposed two alternatives in terms 
of advantageous distortion of the Bromwich contour upon 
which to perform the numerical integration of equation (16). 
The simple parameterized curve has been experimentally found 
to be the more efficient of the two methods and its algo- 
rithm will be developed at this time. The steepest descent 
procedure will be treated in Chapter V. 

The computer program and subprograms which implement 
the simple procedure developed in the preceding chapter are 
all written in the FORTRAN language, using double precision 
arithmetic. They have been tested and debugged on an IBM 
360/67. For machines having a larger mantissa, single pre- 
cision arithmetic might prove to be satisfactory. 

SUBROUTINE VALUE is a user supplied routine which imple- 
ments computation of the g(s) given by equation (11) for the 
desired transform f(s). This subroutine calculates the u(s) 
and v(s) corresponding to a given x and y. 

Either complex or real arithmetic may be utilized in this 
program, although complex arithmetic offers a decided advan- 
tage in convenience and simplicity. However, the user 
should be aware of the potential dangers in unintentional 
exchange of sheets of a Riemann surface when investigating 
a multi-valued function. More will be said about this later. 
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Also, although FORTRAN manuals state or imply that the 
FUNCTIONS DREAL and DIMAG are implemented in all compilers, 
this may not actually be the case. Accordingly it is prudent 
to add the following statements after SUBROUTINE VALUE. 



FUNCTION DREAL (CPLX) 
REAL *8 CPLX ( 2 ) 

DREAL = CPLX(l) 

RETURN 

END 

FUNCTION DIMAG (CPLX) 
REALMS CPLX( 2 ) 

DIMAG = CPLX( 2 ) 

RETURN 

END 



SUBROUTINE CURVE is also a user supplied routine which 
calculates the x and y values for each consecutive incrementa- 
tion of the parameter p of equation (17), thus locating a 
point on the simple parameterized contour. Equation (18) 
provides an illustrative and highly effective example of 
parabolic form which has been employed almost exclusively 
during this investigative work. 

SUBROUTINE INCREP formulates the increasing sequence of 
the parameter p as given by the following equation 

p=p+A (20) 

Limited analysis with other than equally spaced incrementa- 
tion of points along the contour of integration did not ex- 
hibit any enhancing capabilities and was not continued. 

The main program accomplishes the integration of equation 
(.14) using a trapezoidal summation as described by equation 
(16). In addition, the absolute value of the integrand is 
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summed in accordance with equation (19). Also, the number 
of times a sign changes occurs between successive evaluations 
of the addend of equation (16) is recorded. These features 
afford the user an opportunity to monitor potential oscilla- 
tion effects upon the accuracy of the results. Clock time 
and the number of calls to SUBROUTINE VALUE required for 
performing the numerical contour integration are maintained 
as measures of the computational efficiency of the numerical 
inversion procedure. 

Termination of numerical integration occurs when each 
of N successive evaluations of the addend of equation (19) 
has magnitude less than a specified epsilon. We have usually 
used N = 5. This is a very stringent requirement and could 
probably be relaxed resulting in some saving of computational 
time and with but minor loss of accuracy. 

In addition to the termination criterion described in 
the previous paragraph, other program inputs include the 
increment A of equation (18), appropriate numerical values 
of any constants of the function f(s), the starting position 
A of equation (18) and the value of t. 

The program output prints the starting position of the 
integration contour, the values given by equations (16) and 
(19), the final values of x and y at the termination of inte- 
gration, the number of changes of sign between successive 
addends in equation (16), the total number of calls made to 
SUBROUTINE VALUE during the process, and the clock time re- 
quired for numerical integration. 
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As an illustration consider the following example where 



f(s) 



s 



s 2 + 



a 



( 21 ) 



A suitable form for SUBROUTINE VALUE is as follows: 

SUBROUTINE VALUE (X, Y, T, AL, U, V, MMMMM ) 

C THIS IS TABLE ENTRY NUMBER 09 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL *8 ( A , B , D-H , O-Z ) 

IMPLICIT COMPLEX* 16 (C) 

MMMMM = MMMMM + 1 
ZERO=0 . D+0 
CS=DCMPLX ( X , Y ) 

CAL=DCMPLX ( AL , ZERO ) 

CT=DCMPLX(T, ZERO) 

CDEN=CS * * 2+CAL * * 2 
C=CS /CDEN 
CEXP=CDEXP( CS*CT ) 

C=C*CEXP 
U=DREAL(C) 

V=DIMAG( C ) 

RETURN 
END 

If the contour indicated by equation (18) were chosen and 
numerical inversion were performed upon the f(s) given by 
equation (21) for the case where A = 0.3, A =0.1, a = 0.125, 

N = 5 and e = l.D-11, the evaluation of F(2 tt) would yield the 
following output 



AA 


F(2tt) 


| G(2tt ) | X F Y f 


LOSC 


MMMMM 


CO 

o 


7 . 07107D-1 


1 . 32635D+0 -5.46 2.4 


4 
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where AA, LOSC , and MMMMM are the FORTRAN program names for 
the contour starting position, the number of sign changes 
between successive addends, and the total number of calls 
made to SUBROUTINE VALUE, respectively. Similar results 
for the almost 100 functions f(s) tested during the course 
of this research are listed in Appendix B. 



(Note: The parameter a is 

denoted by the FORTRAN 
variable AL. ) 
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IV. AUXILIARY SUBROUTINES 



Chapter III briefly mentioned the potential dangers in- 
herent in the utilization of complex arithmetic within SUB- 
ROUTINE VALUE when the function f(s) is multi-valued. 

If f(s) has branch points or winding points it is impor- 
tant to assure that one remains on the same branch of the 
function as one proceeds along the path of integration. 
Available FORTRAN operations do not assure this and it is 
necessary to employ specially programmed algorithms. The 
present discussion is limited to the two simple cases, rais- 
ing of a complex number to a non-integer power and taking 
the logarithm of a complex number. 

The essence of the matter lies in the FORTRAN operation 
DATAN2(y,x), or its equivalents, which must be relied upon 
to provide the argument of a complex number in polar form. 
This operation always provides a result in the range 
-tt < 0 <_ it, so that if Arg (z) actually passes through the 
values (2n+l)ir, for integer n, the FORTRAN produced value 
of Arg(z) experiences a jump of ± 2 n . If one is determining 
the logarithm of z or a non-integer power of z, shifting of 
branch will take place unless the continuity of Arg(z) is 
restored. 

SUBROUTINES CPOWER and CLOG provide for this continuity 
by monitoring changes in Arg(z) and adding or subtracting 
2nu , whenever it is appropriate to do so, to the value of 
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DATAN2 ( y , x) . As they are presenting written, each provides 
for up to five separate calls, within SUBROUTINE VALUE, to 
either SUBROUTINE CPOWER or SUBROUTINE CLOG. 

These subroutines are listed below, 



SUBROUTINE CPOWER ( ZI N, ZCUT, P, J , NCALL ) 
IMPLICIT REAL *8 (A-H,0-Z) 

COMPLEX * 16 ZIN, ZOUT 

DIMENSI ON T0LDI5) , NCALL ( 5 ) , NPIC5) 

IF(NCALL( J ) . EO. 1) NPI(J)=0 

OH = 0.D+0 

X = DREAL(ZIN) 

Y = DIMAG(ZIN) 

R = DSORT ( (X*X ) + (Y *Y ) ) 

TR0N3 = DAT AN2( Y » X) 

IF(NCALL( J ) .EQ.l) GO TO 5 
PROD = T RONG*TOLD( J ) 

I F ( PROD .GT.OH.OR. X. GT. OH) GO TO 5 
IF (TP ONG.LT .OH) NPI(J) = NPI(J) + 2 
IF(TRONG.G“.OH) N R I ( J ) = NP I ( J ) - 2 

RP s RsJs^rD 
EN = NP I (J ) 

TRITE = TRQNG + ( EN*3.14159265358979324D+0) 
TP = TRI TE * P 
X = R P* DCOS ( T P) 

Y = RP*DSIN(TP) 

ZOUT = DC MP LX( X, Y ) 

T OLD ( J ) = TROMG 

NCALL ( J ) = NCAL L( J ) + 1 

RETURN 

END 



SUBROUTINE CLOG ( ZIN, ZOUT, J, M CALL ) 
IMPLICIT REALMS (A-H,0-Z) 

C OMR LEX* 16 ZIN, ZOUT 

DIMENSION T OL D( 5 ) , MC ALL ( 5 ) , N P I ( 5 ) 

IF(MCALL(J ) .EQ.l ) NPI(J) =0 

OH = O.D+O 

X = DREAL(ZIN) 

Y = DIMAG( ZIN) 

R = DS0RT( { X*X) + ( Y*Y) ) 

TRONG = DA T AN2 ( Y , X ) 

IF( MCALL( J) .EO. 1) GO TO 5 
PROD = TRONG*TOLD(J ) 

I F ( PROD .G T . OH. OR. X. GT.OH) GO TO 5 
IF(TROMG.LT.OH) NPI(J) = NPI(J) + 2 
IFtTRONG. G T .OH) MP I { J ) = NPI(J) - 2 
5 REAL * DLOG(R) 

EN = NP I (J) 

TRITE = +RONG + ( EN*3 .141592653589793240+0) 
ZOUT = DCMPLX(R£AL, TRITE) 

TOLD ( J ) = TRQNG 

MCALL ( J ) = MCALL ( J) + 1 

RETURN 

END 
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If either algorithm is to be utilized by SUBROUTINE VALUE, 
the first requirement is to zero the applicable integer 
array, NCALL or MCALL, within the main program. Then, in 
SUBROUTINE VALUE, the calls to SUBROUTINE CPOWER are of 
the form 

CALL CPOWER ( Z I N , ZOUT , P , J, NCALL ) 

where ZIN is the COMPLEX*16 operand, ZOUT is the C0MPLEX*16 
output, P is the REALMS power, and J is the INTEGER*4 index 
which tells whether this is the first, second,..., call to 
SUBROUTINE CPOWER within SUBROUTINE VALUE. The relationship 
obtained is 

ZOUT = ( ZIN) P (22) 

SUBROUTINE CLOG is used similarly. Calls made from SUBROUTINE 
VALUE are of the form 

CALL CLOG(ZIN, ZOUT, J, MCALL) 

where 



ZOUT = Ln(ZIN) (23) 

J plays the same role as in SUBROUTINE CPOWER. There is no 
P. It should be noted that the array MCALL replaces NCALL. 



24 



V. THE STEEPEST DESCENT CONTOUR 



The theoretical advantages of a steepest descent contour 
have been mentioned oreviously. Consider a contour in the 
complex plane along which v=constant. As shown in Figure 
2, let the £ axis be tangent to the contour with the n axis 
perpendicular to the E, axis in the same sense that y is perpen- 
dicular to x. 







x 



■y REAL A*IS 



Figure 2. Steepest Descent Contour. Contour along 



which v is constant. Axis E, is tangent at 
point P. Axis n is normal. 



Recall equation (11), 



g(s) = e St f(s) = u(s) + iv(s) 



(ID 



By the Cauchy-Riemann relationships 



3u _ _3v . _ _3v _ _3u 

7? " ITn ; 7c ~ ITn 



(24) 



Along a v = constant contour 




( 25 ) 
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Such a v = c = constant contour divides the plane, locally 
at least, into two regions: R^ with v > c and with v < c. 

Consequently along such a contour is one signed, and, 

df| 

therefore, is one signed. u -*■ 0 along the v = c contour, 
in the direction of positive E, . Therefore, there is no 
oscillation in u and, obviously also, there is none in v. 

Since u 0 and v -*■ 0 along a suitable integration contour, 
the selection of c = 0 assures that this results. 

This second program to be discussed does partially follow 
a path of steepest descent. First it integrates vertically 
upward from a chosen point on the real axis (along a Bromwich 
contour) until intersection with a specified v = 0 contour 
occurs. Then it employs a tracking algorithm to obtain suc- 
cessive points on this contour along which it then integrates, 
moving to the left. 

SUBROUTINE VALUE has the same role as that described in 
Chapter III, namely, it calculates the u(s) and v(s) in 
accordance with equation (11) corresponding to a given x and 
y. For the initial part of the algorithm, the integration 
along the Bromwich contour, the value of x is constant, and 
y is increased in increments by Ay. Numerical integration 
in accordance with equation (16) is conducted just as it was 
done in the previous case. 

On the real axis, obviously, v = 0. With an initial 
incrementation in the vertical direction, v assumes a value 
other than zero. The sign of v is recorded and the incrementa- 
tion process is continued until the sign of v changes. This 



26 



constitutes a contour crossing. The procedure continues 
until a predetermined number of such crossings are obtained, 
and then a simple linear interpolation is performed to locate 
the specific ordinate value y for which v = 0. Usually one 
might wish to use the first v = 0 contour above the real 
axis, but the second, third, or other such contour may be 
selected. 

The last (y^,v^) values prior to the desired contour 
crossing and the first (yj> v j) values after crossing are 
utilized in the linear interpolation process to predict YNEW, 
for which v = 0. SUBROUTINE VALUE provides the actual VNEW 
value corresponding to YNEW. Depending upon the sign of 
VNEW, the appropriate coordinate pair (y.,v^) or (y.,v.) is 
replaced by (YNEW, VNEW) and the interpolation is repeated 
with this refinement. Termination of the interpolation occurs 
when YNEW results in a magnitude of VNEW which is less in 
magnitude than a prescribed tolerance limit. 

Once the values of x,y,u, and v have been established 
for the first point on the desired v = 0 contour in the pre- 
vious process, the second, or tracking portion of the 
algorithm commences. 

SUBROUTINE XMARCH is called to provide a new coordinate 
pair (x,y) in the negative sense of the real axis from the 
known, or old, coordinate pair ( x 0 »y o )« This is accomplished 
by the following relationships 

x = x - R cos ( 0 ) (26) 

y = y Q + R sin( 9 ) (27) 
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where (R,9) are the polar coordinates of a search pattern 
measured from the known coordinates (x Q ,y o ) on the v = 0 
contour. The angle 9 is set at zero radians for the first 
call to SUBROUTINE XMARCH. SUBROUTINE VALUE is then called 
to provide the value of v at this point. The sign of v 
establishes upon which side of the v = 0 contour the new point 
lies. 

SUBROUTINE ROTATE is then called to increment the angle 
theta in an appropriate sense. A sweeping search is then 
conducted by SUBROUTINES ROTATE, XMARCH, AND VALUE until 
a sign change is detected for v values between increments, 
thus, indicating a crossing of the v = 0 contour. 

The last (0^,v^) coordinate pair prior to a crossing and 

the first (9.,v.) coordinate pair after a crossing are 
3 9 

utilized in a linear interpolation scheme, analogous to that 
of the first program segment, to refine the location of the 
new point on the v = 0 contour. 

Once the second point on the v = 0 contour is established, 
the tracking algorithm may be enhanced by predicting the 
third point on the contour to lie at the same angle 9 from 
the second point. As before, the sign of v at this new loca- 
tion permits a rotational search to be conducted in the 
appropriate sense. For small values of the radius vector 

R (R=.025 is built into the program, but may be changed by the 
user), this approximation has been found, in general, to be a 
reasonable one. 

The addends of equations (16) and 19) are evaluated at 
each step in the marching process within the main program 
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and, once again, termination of the numerical integration 
occurs when each of N successive addend contributions is 
less, than a prescribed epsilon. 

Required inputs to the main program are the index speci- 
fying the desired v = 0 contour, the starting position on 
the real axis, the increment of the vertical march along the 
Bromwich contour, the termination criterion, any constants 
associated with the function f(s), and the value of t. These 
inputs are represented by the program variables NREG, AA, 
DELY, EPS, (AL,BE,...), and T, respectively. 

The simple parameterized contour offers one decided ad- 
vantage over the steepest descent contour: each call made to 
SUBROUTINE VALUE is productive, in the former case, in the 
sense that it is utilized to compute a contribution of the 
addend of equation (16). The latter algorithm requires calls 
to SUBROUTINE VALUE to perform the linear interpolation pro- 
cedures, which are non-productive in terms of addend computa- 
tion. This represents increased computer time. The present 
investigation has shown that in each of the numerous cases 
which have been examined, the simple parabolic path indicated 
by equation (18) is more efficient and fully as accurate as 
a steepest descent contour. 



29 






VI. SHANKS' ACCELERATOR 



The effectiveness of a family of non-linear sequence to 
sequence transformations in accelerating the convergence of 
(some) slowly convergent sequences and in inducing conver- 
gence in (some) divergent sequences has been reported by 
Shanks [19], If (A n > (n=0, 1 ,2, . . . ) is a sequence of numbers, 

we form a sequence of sequences (A } which, for convenience, 

k , n 

we write as {A(k,n)}. The integer k indicates the "order" 
of the sequence, with A = A(0,n); the integer n indicates 
the position of the term in the sequence. The rule for form- 
ing the sequence (A(k+l,n)} from the sequence (A(k,n)} is 

A ( i , j ) = 0 i f j < 2i 

A ( k+ 1 , n ) = A (k,n)A(k f n-2) - A(k,n-1)_ 2 (28) 

A(k,n)+A(k,n-2)-2A(k,n-l) 

n=2k+3 , 2k+4 , . . . . 

To illustrate the power of this computational device con- 
sider the application of equation (28) and its iterates to 
the first ten terms of the alternating series 

Ln(2) =1 -| + |"| + ••• (29) 

The results are 
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where the tenth partial sum, A(0,10), is correct to only one 
significant figure. However, iterate A(4,10) is correct to 
eight figures since Ln(2) = 0,6931471806. 

In this thesis research, Shanks' accelerator has been 
tested extensively in conjunction with the simple parameter- 
ized contours of integration discussed in Chapters II and 
III. 

Modification of the basic algorithm of Chapter III con- 
sisted of the elimination from the main program of the ter- 
mination criterion for the numerical integration given by 
equation (16), and, instead, performing the integration for 
a finite number of terms. These addends of equation (16) 
were stored in a linear array and, subsequently, transformed, 
within the main program in accordance with equation (28). 

The accelerator did not enhance the rate of convergence 
of the numerical inversion process for any of the cases 
investigated. In fact, if the number of terms of A n , which 
are transformed, becomes sufficiently large so that the 
summation of the addend contributions of equation (16) 
approaches the correct result, the transformation will 
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decrease the accuracy of the result as the denominator of 
equation (28) becomes small. 

Appendix A contains the program listing for the simple 
parameterized contour of integration in conjunction with 
Shanks' accelerator. 

Shanks [19] discusses transformations of higher order. 
These transformations may be worthy of investigation. 



VII. LOCATION OF FUNCTION SINGULARITIES 



The criteria of Chapter II for the distortion of the 
Bromwich contour include the requirement that all of the 
singularities of a function f(s) be to the left of the con- 
tour. Heretofore, the discussion has assumed that the loca- 
tion of these singularities was known. Indeed, this is not, 
in fact, always the case. In this chapter we investigate a 
method of determining the location of the singularities of 
f(s) when such information is either not known or is diffi- 
cult to obtain analytically. 

To see how one can proceed without knowing where the 
singularities are, consider the simple case 



f (s) 



2 2 
s + a 



(30) 



for the particular values of a = 0.125 and t = 2 tt. 

We use a simple parabolic contour given by equation (13) 
with B=1 and with A, which we denote by the symbol x in the 
discussion to emphasize its actual significance, being given 
various values. If we plot the result of the inversion 
versus the starting position x , the result is as shown in 
Figure 3. 

For values of x less than zero, the result is 0.00. 

o 

For x q in the range from about 0.2 to 2.0, the result is 
0.707107. For x q larger than about 2.0, the result is near 
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0.7 but is sensitive to the precise value of x . For x q in 
the range from 0.0 to 0.2, the result exhibits a "spike" 
and a very significant change from 0.0 to 0.707107. 

t.o r 




- 0.5 



-1. -o. v 1. 2. 3- 

AO 

Figure 3. F(t) vs contour starting position for 

cos(at). Value of inverse of s 

for t = 2tt as a function s^+f 125)^ 

of starting abscissa x Q . ' * ' 

We can explain this behavior as follows. For x less 

o 

than zero, the singularities are on the wrong side of the 
contour and we simply get the wrong answer. For x Q in the 
neighborhood of 0.0 to 0.2, the contour passes so close to 
the singularity that the trapezoidal integration scheme is 
not accurate and experiences difficulty in getting an 
accurate value either for the right answer, 0.707107 or for 
the wrong answer 0.0. For x Q in the range from 0.2 to 2.0 
the singularities are on the correct side of the contour and 
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we get the correct answer. For x q larger than 2.0, there are 
numerical difficulties, probably due to oscillation or to 
the effects of the exp(st) factor for s having large real 
part . 

T »Ve have done the same thing with many other transforms 
and have found the same behavior. Our experience leads to 
this suggestion for cases where the user does not know the 
location of the singularities. Perform the calculation for 
a number of different values of x Q and plot the results as 
in Figure 3. The correct result is the ordinate of the 
rightmost plateau, before oscillations or other numerical 
difficulties have had a chance to introduce inaccuracies. 

In the case illustrated in Figure 3, there were only two 
plateaus, one corresponding to all the singularities lying 
on the wrong side of the contour and one corresponding to 
all the singularities lying on the correct side of the con- 
tour. However, with more singularities, one might expect 
to have several plateaus, only the rightmost of which repre- 
sents the correct result. 



VIII. THE HYPERBOLIC TANGENT 



In Chapter II it was stated that great advantage accrues 
from distorting the Bromwich contour as long as the process 
of distortion does not cause the contour to cross over a 
singularity of f(s). In this chapter we investigate a typi- 
cal case in which it is impossible to satisfy this condition. 

The function 



f ( s ) = ^ tanh(^) (31) 

has an infinite sequence of simple poles on the imaginary 

TT 

axis at points s = ±(2n+l)— . By other methods it may be 
shown that this function is the transform of the square wave 
function shown in Figure 4. In this chapter we investigate 
the extent to which our present methodology is capable of 
obtaining this square wave. 
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Figure 4. Square wave function. 

The convergence of the numerical integration has been 
assured by bending the contour to the left in all other cases 



reported in this thesis. Accordingly, we do so in the pres- 
ent instance. This implies that we must pass between two 
poles as the contour bends to the left and all of the in- 
finite number of poles above the point where the contour 
crosses the imaginary axis thus fall on the wrong side of 
the contour and are not included on the left. Thus, in 
effect, our procedure is finding the inverse Laplace trans- 
form of a substitute function which has only a finite number 
of poles coinciding in location with those of the given 
function f(s) and having the same residues. There is some 
reason to hope that the inverse of the substitute function 
will be essentially similar to that of the given function, 
at least for sufficiently small values of t. 

Indeed we do find this to be the case. SUBROUTINE CURVE 
of the algorithm of Chapter II was modified to formulate 
a simple parameterized curve, which utilized the Bromwich 
contour for numerical integration in its initial segment, 
starting from a position s=y on the positive real axis and 
continuing until the elevation of a finite number of poles 
is reached. The second segment of the simple parameterized 
curve which was used consists of the parabola of equation 
(18). This parabolic arc allowed the integration contour 
to cross the imaginary axis between adjacent poles of the 
function f(s). Such a contour is indicated in Figure 5. 

The termination criterion which was used in this case 
is identical to that described in Chapter III. All other 
program input values remained unaltered from the form of their 
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Figure 5. Contour of integration for the hyperbolic 
tangent.. Special contour for inversion of 
f(s) = 75 - tanh(as/ 2 ). It starts at x , 

rises vertically, and then follows a para- 
bolic path between two poles. 



previous presentation with the exception that the integer 
JJJJJ , which selected the number of function poles to be 
encompassed by the contour was added. The program listing 
appears in Appendix A. 

The approximate solutions F(t) obtained by numerical 
inversion of the substitute function for the f(s) given 
by equation (31) are shown graphically below in Figures 6-9. 

Figure 6 is a reproduction of five periods of the square 
wave of Figure 4 with a = 10.0 and the first 100 poles of 
f(s) above the real axis to the left of the contour of in- 
tegration. Figure 7 is a similar picture of the same F(t) 
with a = 5.0 and the first 150 poles of f(s) above the real 
axis encompassed by the contour, with only one period of the 
wave form. Figures 8 and 9 are analogous to Figures 6 and 7, 
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Figure 6. Numerical approximation of five cycles of the 
square wave function with 100 poles to the 
left of the contour. 
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Figure 7. Numerical approximation of one cycle of the 
square wave function with 150 poles to the 
left of the contour. 
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T 

Figure 8. Numerical approximation of five cycles of 
the square wave function with 7 poles to 
the left of the contour. 




T 



Figure 9. Numerical approximation of one cycle of the 
square wave function with 7 poles to the 
left of the contour. 
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respectively; however, only the first seven poles of f(s) 
above the real axis were included on the correct side of the 
contour of integration. 

The reader should be aware that the increment of t used 
in the above series of graphs differed in each case. The 
inaccuracies shown in Figures 6-9 are partially due to the 
assumptions inherent within the employment of our algorithm, 
and, also, partially due to the fact that the plotter simply 
connects points sequentially. This can also be seen in the 
slope of the graphs at t = a, 2a, etc., which is a function 
of the particular At employed and not an inaccuracy of the 
algorithm. 

Efforts to accomplish numerical inversion of the f(s) 
of equation (31) using the steepest descent contour were 
unsuccessful. The difficulty in using this contour appeared 
to be related to an adverse influence from the proximity of 
the function poles to the v = 0 contours near the imaginary 
axis. The algorithm was overtaxed in this vicinity and 
could not track the contour; consequently, further investi- 
gation was not attempted. 

In summary, we conclude that the effect of substituting 
a function which has only a finite number of poles which are 
coincident in location with the first n poles of a function 
f(s), which possesses an infinite number of poles spaced incre- 
mentally along the imaginary axis, produces an inverse 
approximation to the inverse of the given function f(s). 

The accuracy obtained with such a substitution is increased 
as more of the function singularities are encompassed by the 
distorted contour. 
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I X . FACTORS WHICH INFLUENCE ACCURACY 



Heretofore, our discussion has included the treatment of 
contour integration, utilizing a simple parameterized curve 
for that purpose, from two perspectives. The first of these 
methods, discussed in Chapter III, incorporated a termina- 
tion criterion which required a prescribed number, N, of 
successive integrand contributions each to have smaller mag- 
nitude than a specified epsilon. The second of these methods, 
discussed in Chapter VI, performed numerical integration for 
a finite number of points spaced along the contour with sub- 
sequent employment of Shanks' accelerator in an endeavor to 
enhance convergence. 

The first of these methods, we have found, was highly 
successful when employed with the proper point spacing, and 
a very tight termination criterion, namely, a small epsilon 
and large N. The second method has not, in our investiga- 
tions, been found to be successful, and need not be discussed 
further. The issue then becomes one of selecting the proper 
combination of step-size, epsilon, and N to provide an 
accurate and economical result. 

In order to make such an appropriate selection of these 
parameters, it is necessary to examine the factors which in- 
fluence accuracy. These factors include: the step-size along 
the contour, the termination criteria, and the oscillation. 
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Consider a simple case, which we have chosen, along a 
path for which the oscillation is insignificant. The accu- 
racy in such a case is influenced, if the termination 
criteria are sufficiently stringent, only by the step-size 
along the contour of integration. Thus dealing with the 
transform 



f(s) 



~ 2 ~~ Z 

s + a 



(32) 



whose inverse is the cosine function, yields a plot of the 
logarithm of the percent error of F(t) versus the step-size, 

A, as shown in Figure 10. In this case, the contour of in- 
tegration was the simple parameterized curve of equation 
(18) with A = 0.3, B = 1.0, a = 0.125 and t = 2 tt. 

The sequences of points in this plot are all obtained by 
maintaining N constant (N=3) and plotting log^Q percent error 
versus A with epsilon as a parameter. 

Four such parameteric plots appear in the figure below 
corresponding to the values of epsilon (EPS) as shown in 
the legend. 

Clearly the accuracy of the results may be seen to behave, 
in a general sense, in the manner shown in Figure 11, where 

eps 1 < eps 2 < eps 3 < eps 4 . 

When the termination criteria are tight, namely, when EPS 
is small and N is , at least, greater than one, the percent 
error is a narrowly banded, roughly linear function of the 
integration step size, A, over a range from about A = 0.02 
to A = 0.2. This corresponds to EPS^ in Figure 11. 
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Figure 10. Logarithm of percent error versus step size 
for the cosine function. 



If the termination criterion EPS is looser, as evidenced 
in other cases, shown schematically in Figure 11, the accu- 
racy of the result cannot be improved indefinitely by decreas 
ing, A. Rather, there is a point, which differs in each case 
beyond which the accuracy is decreased as the step size is 
decreased. The behavior is complicated by the fact that the 
effects of changes in step-size and changes in termination 
criteria are not themselves disjoint. Making the step-size 
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Figure 11. Schematic representation of accuracy factors. 

Behavior of the logarithm of percent error 
versus step size A as a function of the ter- 
mination criterion EPS. 



smaller results in meeting the termination criteria at an 
earlier point along the path of integration and thus may 
actually reduce overall accuracy. 

One way of dealing with this matter, of course, is to 
revise the termination criteria so as to remove this inter- 
dependence. We simply have not experimented with using 
alternate forms for the termination algorithm. Instead, we 
suggest the following viewpoint to the user who wishes to 
assure that he obtains results with a specified accuracy. 

He should first of all employ a very strict termination 
criterion and use a succession of rather small values of 
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step size, and should also vary the starting point on the 
contour. In this way he should be able to establish an 
evaluation having even greater accuracy than he requires. 
Then, selecting what seems to be a good starting point and 
keeping to it, he should increase the step size until the 
result is no more accurate than he requires. He should 
observe how much saving in computer time (or in the number 
of calls to VALUE) this affords and should not increase step 
size beyond a reasonable point. Then he should loosen the 
termination criterion progressively, measuring the saving 
in computer time versus the possible loss of accuracy, 
stopping short of the point where he cannot rely on obtain- 
ing the required accuracy. 

In many cases where numerical integration is employed 
to obtain a result whose correct value is Y, the numerically 
produced approximation y approximately satisfies a relation 



where a and m are unknown constants. The exponent m may 
frequently be established, once and for all, for the type 
of integration being used, and, for individual integrands, 
the correct value, Y, from two evaluations, y^ obtained with 
A = A-^ and y^ obtained with A = Ag, by using the extrapola- 
tion formula 



y = Y + aA 



m 



( 33 ) 




( 34 ) 
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We have investigated this method of obtaining improved accu- 
racy for our inversions and have found that with the inter- 
dependence between termination criteria and step size, the 
numbers a and m were both rather large and difficult to 
establish - in other words, the extrapolation was not success- 
ful in producing improved values. Our conclusion was that 
it was more efficient from a computational point of view 
simply to use a sufficiently small step size and an appro- 
priately matched set of termination criteria so as to be 
able to obtain an accurate answer without interpolation. 

However, a suggestion for further study of this matter 
is made at the end of the next chapter. 
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X. RESULTS, CONCLUSIONS, AND RECOMMENDATIONS 



This thesis has demonstrated that numerical integration 
along a suitable contour in the complex plane, thus implement- 
ing the complex inversion integral formula, is an effective 
way of obtaining numerical values for the inverse Laplace 
transformation of a given function. 

The success of the procedure , i.e., its efficiency, avoid- 
ance of inaccuracy due to oscillation, and its termination 
after including only a reasonable number of addends in the 
numerical sum, is related to the fact that the contours chosen 
bend to the left as they rise, thus taking advantage of the 
exp(st) factor in the integrand of the inversion formula. 

Although there are some theoretical advantages to what 
we have called the "steepest descent" contour which follows 
a v=0 path, the algorithm which permits following such a path 
requires repeated evaluations of the integrand, many of which 
are not actually used in forming the summands for the numeri- 
cal summation. We have found by investigating almost 100 
cases for which known inverses are available in analytic form, 
that it is more efficient to use a simpler contour than the 
steepest descent path. By using what we have called a simple 
parameterized curve, we have devised a procedure which makes 
use of every evaluation of the integrand. We have found that 
the penalties of not following a steepest descent curve, 
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namely an increase in oscillation of the integrand along the 
path, may be made quite tolerable by using any of a number 
of simple curves. In particular, most of our successful work 
has been with a simple parabola. Our algorithm, however, 
includes the output of information which can serve to alert 
the user to the possibility of inaccuracy due to oscillation 
and thus suggest to him that he might select a different form 
for the path. 

We have even been successful in obtaining the inverse of 
a function which has an infinity of poles spaced along the 
imaginary axis. In doing so we violated the principle 
that all poles of the transform must be to the left of the 
contour. Nevertheless the results are quite satisfactory. 

Our success and accuracy have been so gratifying that 
we venture to suggest that our method may prove to be an 
efficient alternate method for the evalution of exotic func- 
tions for which other methods are slowly convergent or in- 
volve series the coefficients in which are difficult to 
obtain. For example, case 92 in Appendix B shows the suc- 
cessful evaluation of the rather uncommon Struve function. 

The original reason for attempting to employ numerical 
integration in the complex plane as a means of inverting a 
Laplace transformation was that alternate methods employed 
by Hiep and Zargary in conjugated heat transfer problems 
were simply not accurate enough. They led to physically 
impossible results with some temperatures in the media below 
sink temperature and others above source temperature. In 
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this thesis we have not restudied these heat transfer prob- 
lems and recommend that this be done using our methods of 
inversion. There is some reason to believe that the diffi- 
culty encountered by Hiep and Zargary will not be encountered 
if our method is used, but we are not prepared to substan- 
tiate this claim at this time. 

Also it might be of use if the effects of varying the 
available parameters (shape of curve, spacing of points, 
starting point of curves, termination criteria, etc.) were 
to be investigated further so that a user would be better 
prepared to deal with indications of inaccurate inversion 
or inordinate requirements for computer time. Our own 
numerical experiments were invariably so happily successful 
that we have not encountered need for such information. 

For functions f(s) for which it may be difficult to 
locate the singularities, we have shown that varying the 
starting point of a simple parameterized contour permits 
selecting an optimum contour in the sense that one can be 
assured that all singularities are to the left of the 
curve and also that the curve does not reach far enough to 
the right to impose numerical difficulties with large posi- 
tive exponentials. 

At the end of the preceding chapter we indicated that 
we did not find it profitable to employ extrapolation as a 
means of obtaining improved accuracy. However, this is 
probably worth looking at again, and our suggestions for 
doing so are as follows. First employ a termination criterion 
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which is disjoint from the integration step size in the 
sense that the location, along the path of integration, of 
the point where termination occurs, is independent of step 
size. One way of doing this is to employ an epsilon in the 
termination criterion which is itself the result of multiply- 
ing a fixed, input epsilon by the step size. Then the 
extrapolation has a chance of being successful. So as to 
maintain optimum computational efficiency, one should use 
what is called Richardson extrapolation in which the evalua- 
tion points for the larger step size are themselves included 
among those for the smaller step size. 
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APPENDIX A 



LISTING OF COMPUTER PROGRAMS AND SUBROUTINES 

The computer programs and principal subroutines that we 
have prepared and tested during our investigation are all 
listed within this appendix. Section A-l contains the user's 
instructions glossary, flowchart, and program listing for the 
simple parameterized curve. Section A-2 is a similar treat- 
ment: for the steepest descent contour. Section A-3 contains 
the user's instructions and program listing for the simple 
parameterized curve procedure adapted for use with Shanks' 
accelerator. Section A-4 contains the same information for 
adaptation of the simple parameterized curve for use with 
the hyperbolic tangent function described in Chapter VIII. 

It is hoped that sufficient detail has been provided to 
enable a user to adapt one of these programs to his purpose 
in an efficient manner. However, if additional insight is 
required, it may be necessary to refer to the chapter of 
this thesis in which the algorithm is developed. This is 
particularly true in the cases of Sections A-3 and A-4; the 
development in these sections has not been repeated where it 
is equivalent to that of Section A-l. 

Functions DREAL and DIMAG and SUBROUTINES CPOWER and 
CLOG, if they are required, are listed within Chapter III, 
as is a typical example of SUBROUTINE VALUE. Additionally, 
all functions f(s) which have been investigated during this 
research have their applicable SUBROUTINE VALUE listed within 
Appendix B. 



SECTION A-l : THE SIMPLE PARAMETERIZED CURVE 

USER INSTRUCTIONS 

This is the FORTRAN IV program for the numerical inver- 
sion of the Laplace transformation of a function f(s), which 
performs contour integration along a distorted Bromwich con- 
tour in the form of a simple parameterized curve. The fol- 
lowing instructions are provided to assist the user in 
familiarizing himself with the program so that he can adapt 
it to his particular requirements. 

1. SUBROUTINE CURVE is a user supplied group of instructions 
which formulates the numerical integration contour. The 
program listing presently contains a simple parabola,, which 
may be utilized directly if desired. 

2. SUBROUTINE VALUE is a user supplied subroutine which 
calculates the real and imaginary components for g(s) = 
exp(st)*f(s) = u(s) + iv(s). For examples of preparation 
of this subroutine the user is referred to the test cases 
of Appendix B. 

3. The mandatory input variables to the main program are 
as follows: 

(a) AA - REAL*8 starting position (S=AA) on the real 
axis of the contour of integration. (Also denoted 
by A in equation (18) and x q in Chapter VII.) 

(b) AL - REAL*8 value of any constant associated with 

the function f(s). (Note: If the function f(s) 

involves more than one constant, the calling state- 
ments to SUBROUTINE VALUE must be modified accord- 
ingly. 
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(c) DELP-REAL*8 increment of the p-parameter leading to 
spacing of points along the contour of integration. 



(d) T - REAL*8 time for which F(t) is desired. 

4. Input parameters which may be altered by the user, at 
his discretion, are as follows: 

(a) EPS - REAL*8 tolerance parameter for the termina- 
tion of numerical integration. The program listing 
presently contains a value of l.D-11. 

(b) NUMBER - INTEGER*4 number of successive addend con- 
tributions less than the specified value of EPS 
for which numerical integration is terminated. 

5. The output variables which are printed as output by pro- 
gram are defined below: 



(a) AA 


starting position on the real axis 


(b) SUM 


the value of F(t) obtained by numerical 
inversion of f(s). (cf. equation (16)) 


(c) SUMA 


the absolute value of the summation of the 
the addend contributions of the numerical 
integration, (cf. equation (14)) 


(d) X 


the final value of x on the contour at which 
integration was terminated 


(e) Y 


the final value of y on the contour at which 
integration was terminated 


(f) LOSC 


the number of sign changes between successive 
addends encountered during numerical integra- 
tion 


(g) MMMMM 


the total number of calls made to SUBROUTINE 
VALUE during the numerical integration process 
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6. The output of these variables is as follows: 
AA, SUM, SUMA, X, Y, LOSC , MMMMM 



7. The time (seconds) required to perform the numerical inte- 
gration is printed on a line immediately preceding the items 
listed above. 

8. As previously described, the program will perform numeri- 
cal integration along one simple parameterized curve. If the 
user desires, the program may be modified to "sweep" the con- 
tour of integration over a range of starting values, AA, (cf. 
Chapter VII). This may be quickly accomplished in the follow- 
ing manner. 

(a) Replace the value of the integer KSTART from 1 to 
the number of contour integrations desired. 

(b) Replace the REAL*8 value of AA with the following 
statement 

AA = XST ART+ AK *X I NCRM 



where 

XSTART = the smallest starting value of the 
contour of integration the user 
desires 

XINCRM = the Ax between starting positions 

of successive contours of integration 
must also be supplied. 

Thus, the program will perform numerical inversion along 
KSTART contours of integration which differ only in their 
starting positions along the x axis. The user may then 
observe increases or jumps in the value of the output variable 
SUM between successive integrations, enabling him to ascertain 
the locations of function poles. The utility of this proce- 
dure is discussed in depth in Chapter VII. 
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9. If SUBROUTINES CPOWER and/or CLOG are called by SUBROU- 
TINE VALUE in conjunction with use by this program, NCALL 
and/or MCALL must be dimensioned and initialized to zero 
within the main program and added to the calling statements 
to SUBROUTINE VALUE. Additionally, any powers, other than 
integer values, must be entered within the main program and 
passed within the calling statements to SUBROUTINE VALUE. 
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GLOSSARY OF VARIABLE NAMES 



VARIABLE 


DESCRIPTION 


LEGEND 


AA 


starting position on the real axis 
of the simple parameterized curve 


2 


ADD 


computed value of the addend of the 
numerical integration 


1 


ADDA 


absolute value of the addend of the 
numerical integration 


1 


ADDDLD 


saved value of the previous addend 
contribut ion 


1 


AL 


constant associated with f(s) 


2 


DELP 


increment of point spacing along the 
simple parameterized curve 


2 


EPS 


tolerance parameter for termination 
of the numerical integration process 


4 


KOUNT 


integer counter utilized to terminate 
the numerical integration process 


1 


KSTART 


integer counter for the number of times 
the numerical integration is desired 
to be performed, using different start- 
ing positions 


4 


LOSC 


integer counter utilized to record the 
number of sign changes between succes- 
sive addends 


3 


MMMMM 


integer counter to record the total 
number of calls made to SUBROUTINE VALUE 
during the numerical integration 


3 


NN 


integer counter to record the total 
number of addend contributions 


1 


NUMBER 


integer input of the desired number of 
successive addend contributions less 
than a specified epsilon for termination 


4 


P 


parameter of the simple parameterized 
curve 


1 


PROD 


flag utilized to detect oscillation of 
the integrand 


1 
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SUM total addend contribution for the 3 

numerical integration process 

SUMA absolute value of the total addend 3 

contribution for the numerical 
integration process 

T time for which F(t) is to be evaluated 2 

TEST flag for termination of the numerical 1 

integration process 

U real part of the function g(s) 1 

UO previous value of U 1 

V imaginary part of the function g(s) 1 

VO previous value of V 1 

X real part of s 1 

XO previous value of X 1 

Y imaginary part of s 1 

YO previous value of Y 1 

LEGEND 

1. No action required by the user 

2. Mandatory user input 

3. Appears as program output 

4. May be altered by the user at his discretion 
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non oo<"> ooooo ooo° oonoooo 



C 

C 



THIS IS THE PROGRAM TC PERFORM NUMERICAL INVERSION CF THE LAPLACE 
TRANSFORMATION OF A FUNCTION USING THE SIMPLE PARA METERI Z ED CURVE 
IMPLICIT REALMS (A-H,C-ZJ i rprM,^ 

CALL ERRSET ( 20 7, 25 6, - 1 » 1, 1 , 2 OS ) LLCL.'Jl 

NUMBER = 5 £ 

KSTART = 1 ^ 

TWOPI = 4. D+Q* DA R SIN ( 1*0+0) ^ 

AL = 1.250-1 ? 

4 t%*n*4mkH4k**44**4U**n* *************************** *4* ******* ****** 
* * 

« THE OUTER lD0P MAY 6E UTILIZED TC SWEEP THE CCNTCUR STARTING * 
^POSITION OVER A DESIRED RANGE AS DESCRIBED IN THE USER INST. * 

* a 

\ 

DO 400 KKK=1* KSTART 

ADO = 0.0+0 

ADDA = C.D+O 

SUM = 0 .00 

SUMA = C.D+O 

LOSC = 0 

MMMMM = 0 

AK = KKK ^ 

AA = 3.D-1 

T = TWOPI *. ^ 

P = O.D+O - 

NN = 0 
KCUNT =0 

EPS = 1.0- H £ 

OELP = 1 .0-1 ^ 

CALL CURVE (AA,P f XO,YCJ ^ 

CALL VALUE <X3 » Y3 , T , A L, U0» VO, MMMMM) 

TIME = O.D+O 
CALL SETIME 

3 $ $$$$$$*$$$$*>$ $*£$$$ 

$ $ 
S THE FOLLOWING LOOP PERFORMS THE NUMERICAL INTEGRATION S 

200 NN = NN+1 

CALL INCFEP (P,OELPJ 
CALL CURVE (AA ,P.X, t ) 

CALL VALUE <X,Y , f ,AL ,U,V , MMMMM } 

ADO » ( ( V+VO) *1 X- XO) + ( U+UOJ * (Y-YQJ ) /TWQPI 

IF (NN • EO. 1 ) ADDOLD^AOO 

PROO * ACG-A003LD 

IF (PRCO.LT.G.j+OJ l0SC=L0SC+1 

ACDA = CABS ( ADD J 

SUM = SUM+ADO 

SUMA = SUMA+AODA 

* * 

* THIS IS THE TERMINATION CRITERION SEGMENT OF THE ALGORITHM * 

♦ ■ * 

TEST = DA3S (ADD) 

IF (TEST .LT .EPS ) KOUNT=KCUNT+ 1 
IF (TEST.GE.fcPSJ KOUNT =0 
CALL GET I u E (1ST) 

EL * OFLGAT (IET J * 2 . o C- 5 
TIME = T I M E + El 
CALL SETIME 

IF (KOUNT. EO. NUMBER) GO TO 300 

* * 

XO = X 
YO = Y 
UC = U 
VC = V 

ADDOLD = ADD 
GO TO 2 00 

s s 

**S3$SSl*$$S$$$$S$S*SS*SSS3$$$3$S*S$SS*SSJS$S3SI$$$iSSS$$i$$S5S$$$ 
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170 
180 
190 
200 
210 
22 C 
230 
240 
250 
260 
270 
260 
290 
300 
210 
320 
33 0 
°4Q 
-sO 
36 0 
370 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
49 0 
500 
510 
520 
53 0 
540 
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560 
570 
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600 
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62 0 
630 
640 
6 50 
66 0 
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680 
690 
700 
710 
720 
730 
740 
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OUOo 



300 CONTINUE 

WRITE <6,5001 NN,TIMt 

WRITE (6,600) Aa, SUM # SUMA,X,Y,L0$C, MMMMM 
400 CCNTINUE 

* # 
STOP 

c 500 FCRMAT ( IX, 1 TOTAL ELAPSEC TIME AFTER • , 15 , * ITERATIONS IS *,620.5) 
600 FORMAT ( IX, IP 5E 20 . 5 , 2 I 8 , / ) 

END 

SIBROUT INE CURVE (A«P«X«Y) 

IMPLICIT REALMS ( A-h,Q-Z ) 

X = A- ( P*P) 

Y = P 

RETURN 

END 

SIBROUT INE INCREP (P,OELP) 

IMPLICIT REAL*d (A-H.O-Z) 

P = P+QELP 

RETURN 

END 



750 
76 0 
770 
780 
790 
800 
810 
820 
830 
840 
850 
66C 
870 
880 
10 
20 
30 
40 
50 
60 
10 
20 
30 
40 
50 
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SECTION A-2 : THE STEEPEST DESCENT CONTOUR 

USER INSTRUCTIONS 

This is the FORTRAN IV program for the numerical inver- 
sion of the Laplace transformation of a function f(s), which 
performs contour integration along the steepest descent path 
in the complex plane. The following instructions are pro- 
vided to assist the user in familiarizing himself with the 
program so that he can adapt it to his particular require- 
ments . 

1. SUBROUTINE XMARCH is a subroutine which tracks a speci- 
fied v=0 contour (steepest descent) in the complex plane, 
once this contour has been intersected for the first time. 

As shown in the program listing, the following two equations 

X = XHOLD - l.D-l*DCOS( THETA) 

Y = YHOLD + l.D-l*DS IN (THETA) 

march the numerical integrations in the favorable direction 
along the contour. The equations are parametric in THETA, 
which is the search angle provided by SUBROUTINE ROTATE. 

The l.D-1 is a polar radius which may be altered to a larger 
or smaller increment by the user, or accepted as it appears. 

2. SUBROUTINE ROTATE increments the angular search angle 
THETA in order to locate a new point on the v=0 contour. The 
angle is presently incremented by 5.D-2 radians with each 
call to the subroutine. This may be altered as required. 

3. SUBROUTINE VALUE is a user supplied subroutine which cal- 
culates the real and imaginary components for g(s) = exp(st) 

x f(s) = u(s) + iv(s). For examples of preparation of this 



65 






subroutine, the user is referred to Appendix B. 

4. The mandatory input variables to the main program are 
listed as follows: 

(a) AA - REAL*8 starting position (s=AA) on the real axis 

(b) AL - REAL*8 value of any constant associated with 

the function f(s). (Note: If the function f(s) 
involves more than one constant, the calling 
statements to SUBROUTINE VALUE must be modified 
accordingly. ) 

(c) DELY - REAL*8 increment of the point spacing along 

the Bromwich contour. (Note: to alter the point 
spacing along the v=0 contour, the user is 
referred to 1 above.) 

(d) NREG - INTEGER*4 input value which specifies the 

first, second, third, etc. v=0 contour above 
the positive real axis, which is to be followed 
in the numerical integration process. 

(e) T - REAL*8 time for which F(t) is to be evaluated. 

5. Input parameters which may be altered by the user, at 
his discretion, are listed below as follows: 

(a) EPS - REAL*8 tolerance parameter for the termina- 

tion of numerical integration. The program 
listing presently contains a value of l.D-11. 

(b) NUMBER - INTEGER*4 number of successive addend 

contributions less than the specified value of 
EPS for which numerical integration is terminated. 
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6. The output variables which are received from the program 
are defined below: 



(a) MMMMM - the number of calls made to SUBROUTINE VALUE 

(b) SUM - the total addend contribution for the numeri- 

cal integration 

(c) SUMA -the absolute value of the total integrand 

contribution for the numerical integration 

(d) EL -the computation time in seconds required for 

the process 

7. A sample program output is as follows 



694 

0. 724724D+0 0.324972D+1 

0 . 707107D+0 0 . 326734D+1 

4027 

0 . 71090D+1 



MMMMM 
SUM, SUMA 
SUM , SUMA 
MMMMM 
EL 



at end of vertical 
section 
at end of 
computation 
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GLOSSARY OF VARIABLE NAMES 



VARIABLE 


DESCRIPTION ] 


LEGEND 


AA 


starting position on the real axis of 
the Bromwich contour segment of the 
integration 


2 


ADD 


computed value of the addend of numeri- 
cal integration 


1 


ADDA 


absolute value of the addend of numeri- 
cal integration 


1 


AL 


constant associated with f(s) 


2 


CHECK 


flag to ascertain when a v=0 contour 
crossing has occurred 


1 


DELY 


increment of point spacing along the 
Bromwich contour segment of numerical 
integration 


2 


EPS 


tolerance parameter for termination of 
numerical integration 


4 


EXAM 


flag to ascertain when a contour crossing 
has occurred 


1 


J 


integer counter used internally for test- 
ing contour regions 


1 


KOUNT 


integer counter utilized to terminate 
the numerical integration process 


1 


KSIGN 


integer counter utilized to record the 
number o f v=0 contour levels crossed 


1 


LEVEL 


integer flag assigned to regions of 
positive and negative v 


1 


MMMMM 


integer counter for the number of calls 
made to SUBROUTINE VALUE during each 
segment of the numerical integration 
process 


3 


NREG 


the desired v=0 contour level above the 
positive real axis which is to be tracked 


2 


NUMBER 


integer input for the number of succes- 
sive addend contributions less than the 
tolerance parameter required for termi- 
nation 


4 
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SUM 

SUMA 

T 

TERM 

TEST 

THETA 

THETAO 

THNEW 

THONE 

THTWO 

U 

UHOLD 

V 

VHOLD 

VNEW 

VOLD 

X 

XHOLD 

Y 

YNEW 



total addend contribution for the 3 

numerical integration process 

absolute value of the total addend 3 

contribution over the contour of 
integration 

time for which F(t) is to be evaluated 2 

termination parameter for the absolute 1 

value of the addend 

Flag to ascertain when a contour 1 

crossing has occurred 

search angle used in the contour 1 

tracking algorithm 

previous value of THETA 1 

predicted angular location of a new 1 

point on the contour from a known 
previous point • 

coordinate utilized in linear inter- 1 

polat ion 

coordinate utilized in linear inter- 1 

polat ion 

real part of g(s) 1 

previous value of U 1 

imaginary part of g(s) 1 

previous value of V 1 

actual value of the imaginary part of 1 

g(s) for the predicted v=0 location 

previous value of V 1 

real part of s 1 

previous value of X 1 

imaginary part of s 1 

predicted elevation of the desired v=0 1 

contour level 
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LEGEND 



1. No action required by the user 

2. Mandatory user input 

3. Appears as program output 

4. May be altered by the user at his discretion 
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71 




72 





73 







RETAIN THE 
X, Y, U, AND V 
VALUES OF 
THE POINT 
ON THE COUNTER 



± 

COMMENCE THE 
TRACKING ROUTINE 
TO LOCATE NEW 
POINTS ON THE 
COUNTER 



V 



1100 
CALL 
X MARCH 
ROTATE 
AND 
VALUE 







-4 

v 



ESTABLISH 
THE SIGN 
OF V 

SET OLD V = V 



74 





75 



I 





76 




1000 

A 

LT. 5 
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c 

c 



TRANSFORM OF A 



400 

5 00 
600 



700 

800 

900 



INTEGRATES ALONG ThE 
THE UESIREO V= 0 
REFINE TFE LOCATION. 

LEoEND 

^ 2 

^ o 



> 



200 



300 



THIS PROGRAM COMPUTES THE INVERSE LAPLACE 
FUNCTION USING Trtt STEEPEST DESCENT CONTOUR 
IMPLICIT RE AL*8 (A-H,0-Z) 

CALL SETIME 

CALL ERRSET < 207 , 256 , - 1 , 1 , 1 ,2 0 9 ) 

**$**£*$#$$ jfc***:*:**#**#*****:**::*************#**#****#*#******* 
* 

* THIS SEGMENT OF THE PROGRAM NUNERICALLY 
* BROMW I CH contour until it INTERSECTS WITH 
^CONTOUR. LINEAR INTERPOLATION IS LiSEO TO 

* 

TWOPI = 4.D+0+0ARSINI 1 .C+0) 

T = TWOPI 

KKKKK = 0 

AL = 1 .250-1 

J = 1 
MPMMM = 0 
SLM = 0.0+0 
SUM A = 0.0+0 

A A = 3.C-1 

X = AA 

Y = 0.0+0 

Y POLO = 0.0+0 

OELY = 2.50-3 

NREG = 3 

KSIGN = 0 

LEVEL = -1 

CALL VALUE ( X , Y , T , AL , LHOLO , VHOLO , RMMMM) 

Y = Y+OELY 

CALL VALUE I X , Y , T , A L. L , \ ) 

J = J ♦ I 

IF (J.EC.2) SIGN = V/CAESIV) 

IF (J.GT.2) TEST = V«SIGN 
IF (TEST .lE.O .0+0 I GO Tu 400 
ACO = I (U+UHOLOJ*<Y-YhCLC) J/TWQPI 
SLM = SUM+AOO 
ACOA = CA8S < ADO I 
SUMA = SUMA+AOOA 

IF (KKKKK. EC.1) WRITE (6,1800) AOO.AOOA 
UHOLO = U 
VHOLO = V 
YHOLO = Y 
GC TO 200 
KSIGN = KSIGN+1 
LEVEL = -LcVEc 

IF (KSIGN. cC.NREuJ GO TC 500 
SIGN = -SIGN 
GC TO 300 






* 



YCNfc = Y-QELY 
Y7W3 = Y 
VCNE = VHQLC 
VTWO = V 

YNE W = ( VQNE* ( YTWQ- YGNE) / ( VONE- VT WC ) ) ♦ YCNE 
CALL VALUE (X ,Y \I£W,T, A L » UN E a , VN E W , MMMMM ) 

IF COABS<VNtWJ. LT.1.0-1U GO TO 9C0 

TEST = VNEW*SIGN 

IF (TEST. GT .0.3+0 J GO TO 700 

Y TWO = YNEW 

VTWO = VNEW 

GO TO 800 

YQNE = YNEW 

VONt = VNEW 

CONTINUE 

GC TO 6 C C 

AGO = C C UN EW + UHOLQ) *( YNEW- YHQLO ) ) /TWOPI 
ACCA = CAdS ( A JO i 
SLM = Sl.M+AOO 
SUMA a SUMA+AOOA 

if f igKfKK-fC.l ) WRITE (6,1800) AOC,AOOA 
WRI TE ( 6,2100) MM MM M 
WRITE (6,2000) SUM, SUMA 

******* *****************4***********#**********^^^^^ *** * ***** 



10 
2C 
30 
40 
50 
60 
7 C 
80 
90 
IOC 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
26 0 
27C 
280 
290 
300 
210 
320 
33 C 
340 
350 
260 
370 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
52 0 
530 
540 
550 
560 
570 
580 
590 
600 
610 
620 
63 0 
640 
650 
66 0 
670 
680 
690 
700 
710 
720 
73 C 
740 
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1000 

1100 

1200 

1300 

1400 

1500 

1600 



1700 



$*i $$SS UJS $$$$$$$ $$!$$$$$ $$$$$$$ $$$$$$$$$$$$$ 

$ s 

$ THIS SEGMENT OF THE PROGRAM PERFORMS NUMERICAL INTEGRATION S 



SALQNG THE V = 0 CONTOUR. 
$ 

K0UNT = 0 
NUMBER = 5- 
EPS = 1.0-11- 



LEcti\C 



THETA = C.0+0 

V = YNE Vi 
J - 1 
XHQL 0 = X 
YhOLO = Y 
UHOLO = UNEW 

V POL 0 = VNEW 
MMMMM = 0 

J = J+ 1 

CALL XM ARC H ( XHOL 0 , YHOL D , THE T A , X , Y) 

CALL VALUE ( X ,Y , T , A L , U , V , MMMMM ) 

IF (J.EG.2) SIGN = V /DAS S ( V ) 

IF (J .GT .2 ) EXAM = V-SIGN 
IF (EXAM.LE.0.0+01 GO TC 1200 
VOLO = V 
ThETAO = THETA 

CALL ROTATE < Lt VEL , SI GN , THETA J 

GO TD 11C0 

ThCNE = THETAO 

THTWO = THETA 

VONE = VOLO 

VTWQ = V 

THNEW = <UTHON£-THTWCJ/CVTWO-VCNE))*VONE)*THCNE 
CALL XMARCH (XHOLOt YHOLO, THNEW, X,Y) 

CALL VALUE (X , Y ,T ,AL , UN EW,VNEW , MMMMM J 
IF (DABS! VNEWJ.LT.1.0-1U GO 1600 
CHECK = VN EVASION 
IF ( CHECK . GT .0 .0+0 J GO TO 1400 
THTWO = THNEW 

V TWO = VNEW 
GO TO 1500 
TkONE = THNEW 

V CNE = VNEW 
CCNTI NUE 

GO TO 1300 

ACQ = { (UN EW+UHJL Di *< Y-YhOLD) + (VNEW+VHOLO I *(X-XHOLD) ) /TWOP I 
AGO A = DABS ( ADO ) 

SUM = SUM + A JO 
SUMA = SUM A+AOOA 

IF (KKKKK.EQ.1J WRITE (6,1800) ADC,AOOA 

***************************** ****4X**iX* 
4 4 

4 THIS IS THE TERMINATION CRITERION SEGMENT OF THE ALGORITHM. * 

4 4 

TERM = DABS(AOO) 

IF (TERM. LT. E°S J KOUNT = KQUNT +1 
IF (TERM.GE.EPS1 K Q U N T = 0 
IF (KOUNT. EO. NUMBER) GO TO 1700 

4 4 
444444444444*4444444444444444444**4444***44*4444444*44444*444444*4 

THETA = THNEW 
EXAM = -EXAM 
GO TO 1000 

WRITE (6,20001 SUM, SUMA 
WRITE ( 6 , 2 1 CO J MMMMM 
CALL GET I M E ( I ET ) 

EL = DFLGATUcT J^2.6C-5 
WRITE (6,19001 EL 

5 S 

STOP 



I 

75 0 

76 0 
770 
780 
790 
800 
810 
82 0 
830 
840 
85 0 
860 
87 0 
880 
890 
900 
910 
920 
930 
940 
950 
960 
970 
980 
590 

1C00 

1010 

1020 

1030 

1C40 

1050 

1C60 

1070 

1080 

ir?o 

1100 

1110 

1120 

1130 

1140 

1150 

1160 

1170 

1180 

1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1350 

1400 

1410 

1420 

1430 

1440 

1450 

1460 

1470 

1480 
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c 

1800 

1900 

2000 

2100 



100 
2 00 



2 . 

4. 

5. 



( 1X,2E30. 10,/) 

( IX, •elapsed time 

(IX ,2E20.o) 

( IX, Id) 



FORMAT 
FORMAT 
FORMAT 
FORMAT 
END 

SUBROUTINE XMARCH ( XHCLC , YHOLO , TH ET A ,X , Y ) 
IMPLICIT REALMS (A-H,0-Z) 

X = XHOLD-2 .5D-3*CC3S (Tl-ETA) 

Y = YH0LQ+2.50-3*DSIN< THETA) 

RETURN 

END 

SUBROUTINE FOTATE (LEVEL ? SIGN, THETA ) 
IMPLICIT 0 E AL* 8 (A-rl,0-Zf 
IF (LEVEL. NE-U GO TO ICO 



DURING THE COMPUTATION IS f ,E20.5J 

LqGLMU 



5 

5 



I F < SIGN.UT.G.O + Gi 
IF (SIGN. GT. 0.3+0) 
GC TO 200 

IF ( SIGN. LT. 0.0+0) 
IF (S IGN.GT .0 .0+0 ) 
RETURN 
END 



THETA=THETA-5.C-2 

THETA=THETA+5.D-2 

THETA=THETA+5. C-2 
THETA-THETA-5.C-2 



1490 

1500 

1510 

1520 

1530 

1540 

10 

20 

30 

40 

50 

60 

10 

20 

30 

40 

50 

60 

70 

80 

90 

100 



LEGEND 



Mandatory user output. 

May be altered by the user at his discretion. 

Polar radius R = .025 built into program, may 
be changed by user. 
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SECTION A-3 . SHANKS' ACCELERATOR 
USER INSTRUCTIONS 

The algorithm of this section is very similar to that 
of Section A-l. Only those features which differ have been 
discussed in this instruction. 

1. Contour integration is performed using a simple para- 
meterized curve for a finite number of points on this con- 
tour. This finite number of terms must be prescribed by 
the user. The integer variable M is used for this purpose. 

2. The dimension statements at the beginning of the program 
must provide storage space in linear arrays which are equal 
to or larger than the value assigned to M. (These quanti- 
ties will subsequently be used with Shanks' accelerator. 

Note particularly that the array is doubly indexed.) 

3. The integer variable N specifies the number of iterates 
of Shanks' accelerator which are desired by the user. How- 
ever, the assigned value of N cannot be greater than 6 as 
the program is now written. 

4. A sample program output appears in Chapter VI of this 
thesis. Column 1 represents the natural result obtained 
after M evaluations of the addend, without application of 
the accelerator. 
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noon on oo no oo onoooooo oonooooo opo 



THIS PROGRAM CALCULATES THE INVERSE LAPLACE TRANSFORM OF A 
FUNCTION USING A FINITE NUMBER OF TERMS ALONG A SIMPLE 
PARAME TERI ZED CONTOUR CF INTEGRATION WITH S HANKS 9 ACCELERATOR 
IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION A ( 10 , 1 JO J • ACO(IOO) 

0 1 MENS 1 CN XX 100) ,YY( 100) ,UU( 100) ,VV(10 0) 

CALL ERRSET ( 20 7, 25o, - 1, 1, 1, 209) 

MMMMM = 0 

TkOPI = 4.0+0*DARSIN< 1.D+0) 

AL = 1.250-1 
T = TWOPI 
AA = 3.D-1 
P = 0*0+0 
M = 20 
NN = 0 
N = 5 

CALL CURVE <AA,P,X0,Y0) 

CALL VALUE ( XO, YJ , T ,AL , LQ , VO , MMMMM ) 

SSSSSSS $$$$ $$$$$$ SSSSii $ $$$$$$$$$$$*$$$ $$$$$$$$$$$$5$$ 



$ s 

$ THE FOLLOWING SEGMENT OF THE PROGRAM PERFORMS THE CONTOUR $ 

$ INTE GR ATI ON FOR A FINITE NUMBER CF TERMS USING THE SIMPLE S 

^PARAMETERIZED CURVE. S 

$ $ 



200 NN = N.N + 1 

CALL INCFEP (P) 

CALL CURVE (AA,P,X,Y) 

CALL VALUE ( X * Y * T * A L » U * V * M MM M M ) 

ACD(NN) = < (V+VQ) XX-XO J + (J+JO )*< Y-YO))/TWOPI 

XO = X 

YO = Y 

UC = U 

VC = V 

XX(NN) = X 

YY(NN) = Y 

UU( NN J = U 

VV(NN) = V 

IF (NN .EO.M ) GO TO 300 
GC TO 200 
300 CONTINUE 

$ $ 

******* *** * ************ ************ ********** ***** **************** 

* * 

♦ THIS SEGMENT OF THE PROGRAM IMPLEMENTS SHANKS 1 ACCELERATOR IN * 

♦ ACCORDANCE WITH EQUATION 28 OF CHAPTER 6 OF THIS THESIS* * 

♦ * 
B = 1*0+0 

EFS = 1 .0-12 



DC 400 1=1, N 



DC 400 J=L,M 
400 A ( I , J ) = 0.0+0 



1 = 1 



500 



00 500 J=1 , M 
JM =’ J-l 
8 = ADD ( J ) 

IF (J.EC.I) 
IF (J.GT.l ) 
CONTINUE 



A ( 1 ,1) 

A( 1, J ) 



e 

B+A( 1 , JM) 



I 

10 
20 
30 
40 
50 
60 
70 
80 
90 
100 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
270 
280 
290 
300 
310 
320 
330 
340 
35 . 
360 
370 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
520 
530 
540 
550 
560 
570 
580 
590 
600 
610 
620 
63 0 
64 0 
650 
660 
670 
680 
690 
700 
710 
720 

73 C 

74 0 
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noon oooo oo 



DO 600 1=2 »N 
IM = 1-1 
K = 2*1-1 



600 



DO 600 J=K, M 
JH = J-l 
JMM = J-2 
A(1,J) = 1 « D+0 

DEN = A(IM # J)+A(I M, JMM)-2.C+0*A(I^ f JM) 

IF (DABS (DEN) .GT.EPS) A ( I » J) = ( A ( I M , Ji *A ( I ^ , JM M -A ( I P . JM )**2 ) / 0 EN 
CONTINUE 



DC 700 J=2,M 

700 WRITE ( 6t 80 C ) J , ( A ( I , J ) , I =1 f N ) 

* * 
******* 4 ********$** ****** ********************** 4 * ************* 4**4 

STOP 



800 FORMAT (5X , 13 , 6 FI 8 .14 i 
END 

SUBROUTINE CURVE (A,P,X,Y) 
IMPLICIT REAl*8 (A-H.C-2J 
X = A- ( P *P ) 

Y = P 

RETURN 

ENO 

SUBROUTINE INCREP (P) 

I KPLICI T REALMS I A-H » C-Z i 

P = P + 1.0-1 

RETURN 

ENO 



750 
76 0 
770 
780 
790 
800 
810 
820 
830 
840 
850 
860 
870 
880 
890 
900 
910 
920 
93 0 
940 
950 
960 
970 
980 
990 
1C00 
1010 
10 
20 
30 
40 
50 
60 
10 
20 
30 
40 
50 
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SECTION A-4. SPECIAL CONTOUR ALGORITHM FOR CASES AS 
TREATED IN CHAPTER EIGHT 

USER INSTRUCTIONS 

The integration contour used in this section is a vertical 
path to a specified value of y, followed by a simple parabola. 
The algorithm of this section is very similar to that of 
Section A-l. Only those items which differ are discussed in 
this section. 

1. Contour integration is performed in the initial segment 
of the program along the Bromwich contour. The integer vari- 
able JJJJJ specifies the number of poles of the function 

1 ct s 

f(s) = — tanh(-?;-) above the real axis that are to be enclosed 
v ' s 2 

to the left of the distorted Bromwich contour, which is used 
in the second segment of the program. 

2. The variable HALT computes the elevation above the real 
axis corresponding to the input value of JJJJJ for the 
specified hyperbolic tangent. 

3. The variable FLAG is internally assigned a value of zero 
or one, corresponding to the first or second segments of the 
contour of integration which the algorithm employs. No 
action is required by the user. 

4. The integer variable KSTART is utilized in the outer loop 
of the program to increment the value of t for which F(t) 

is computed. 

5. The linear arrays STORE and STASH are used to convert the 
values of F(t) and t to single precision, which is required 
for the graphics package associated with the IBM 360/67. 
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6. SUBROUTINE VALUE contains the real arithmetic coding of 
the function g(s) for the hyperbolic tangent function of 
Chapter VIII of this thesis. Standard trigonometric identi- 
ties for the hyperbolic functions of a complex argument have 
been utilized. 

Note: This procedure is obviously keyed to the particular 

function f(s) = — tanh(-^-). The user must modify it for 
other functions. One obvious modification which would 
generalize the procedure would simply be to specify rather 
than calculate the value of HALT, i.e., input HALT rather 
than «J J J J . 
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noon 



THIS PROGRAM COMPUTES THE INVERSE LAPLACE TRANSFORM OF A 
FUNCTION USINl, THE SIMPLE PA R AME TE R I Z EO CONTOUR OF INTEGRATION 
WITH MODI F I CAT I ON TO VARY THE tLcVATION FROM WHICH THE CURVE IS 
LAUNCHED* THE BROMWICH CONTOUR IS USED IN THE INITIAL SEGMENT. 
ILLICIT REALMS (A-h,0-Z) 

REAL*4 STORE, STASH 
DIMENSION STORE! 50), STASH! 50) 

CALL ERRSET (ZO 7 , 25 6 * - i , 1, 1,209) 

TWOPI = 4.D+0*0ARSi,N( 1 .0+0) 

AL = 1 .D + l 
JJJJJ = 15 
KKKKK = 0 
NUMBER = 5 
KST ART = 50 

00 400 KKK = 1, K START 
FLAG = O.D+O 

HALT = !DFLOAT( JJJJJ) *TWQPI)/ (2*D+0*ALJ 
ACD = 0 .0+0 
ACDA = 0 .0+0 
SLM = 0.0+0 
SUMA = 0 .0+0 
LCSC = 0 
MMMMM = 0 
AK = KKK 
AA = 3.0-1 
T = 5.D- l+AK*l.D+0 
P = 0.0+0 
KCUNT = 0 
NN = 0 

EPS = 1 .0-1 1 
0 EL P = 1.0-1 

CALL CURVE !AA,P, HALT, FLAG, XO,YO) 

CALL VALUE <X0 ,YQ ,T , A L , UC , V 0, MMMMM) 

TIME = C.0+0 
CALL SETIME 
200 NN = NN+1 

CALL INCREP ( HA LT , OE LP , F ,F L AG ) 

CALL CURVE (AA, P, HALT , FLAG, X, Y ) 

CALL VALUE (X,Y,7,AL,U,'v, MMMM M ) 

ADO = < ( V+VO) *( X-XJ )+( U + UO) *( Y- VO ) > /TWOPI 

IF (NN.EG.l ) AO JOLD= AOC 
PROD = AD0*A0D0LJ 
IF (PROD.LT .0.0 + 0) L0SC=L0SC+1 
ADDA = DABS (ADD) 

SLM = SUM+ADO 
SUMA = SUMA+AOOA 
TEST = D ABS ! ADD) 

IF ( TEST.LT.EPS) 

IF (TEST .GE .EPS ) 

CALL GETIME !iET) 

EL = DFLOATd £T)*2. 60- 
TIME = TIME+EL 
CALL SETIME 

IF(KOUNT .EQ. NUMBER) GO 
XC = X 
YO = Y 

uo = u 

VC = V 

AOOOLO = AO 0 
GO TO 200 
CCNT INUE 
WRITE (6,500) 

WRITE (6,600) 

STORE ( KKK ) = 

STASH (KKK) = 

CONTINUE 

CALL PLCTG(ST0RE,STASF,50,i, 1,0, 

16 • Oi 6. 0) 

CALL PLOT (0 .0,0 .0, 999 ) 

STOP 

500 FORMAT (IX, ’TOTAL ELAPSED TIME AFTER M5,» ITERATIONS IS *,£20.5) 



300 



400 



KOUN 7 = K CUNT+ 1 
KJUNT =0 



70 300 



NN , T I M£ 
T, SUM, 
SNGL (T ) 
S.NGL (SUM) 



SUMA, X, Y, LOSC, MMMMM 



f T • • 1, *F(T) *,4,0.0, 55. 0, -2.0, 2. C, 
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30 
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60 
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80 

90 



100 

120 
130 
140 
150 
160 
170 
180 
190 
20G 
210 
220 
23 0 
240 
250 
260 
280 
290 
200 
320 
330 
340 
350 
360 
270 
380 
390 
400 
410 
420 
430 
440 
450 
460 
470 
490 
500 
510 
520 

540 

550 

560 

570 

580 

590 

600 

610 

63 0 
640 
650 



660 

670 

680 

690 
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600 FORMAT ( IX , IP 5E20 .5* 2 I 8t / ) 

E NO 

SUBROUTINE CURVE ( A A , P , H ALT ,F L AG , X , Y) 
IMPLICIT REALMS (A-h,0-Z) 

IF IP. LT.HALT.ANO. FLAG. NE. 1.0*0) CO TO 
X = AA-UP-hALT )MP-rtALT) ) 

Y = P 

GC TO 200 
100 X = AA 

Y = P 

2 00 RETURN 
END 

SUBROUT INE INCREP ( HALT, DEL P, P ,F LAG ) 
IMPLICIT REAL*3 <A-H,C-Z) 

P = P + DEL 0 

IF (P.EQ.HALT) FLAG=1.D*0 

RETURN 

END 

SUBROUTINE VALUE ( X , Y » T , AL» U, V,MMMMM) 
IMPLICIT REAL*3 IA-H,C-Z) 

MMMMM = MMM MM* 1 
T EMI = <AL*X) 

TEM2 = ( AL*Y ) 

TEM3 = DSINH(TEMI) 

TEM4 = CSIMTEM2) 

TEM5 = DCCSHITcMl)*3CCS(TEM2) 

TEM6 = TEM3/TEM5 
TEM7 = TEM4/TEM5 
T EM 8 = (TEM6*X)*( TEM7*Y) 

TEM9 = (TEM7*Xj-< TEM6 + Y ) 

T EMI 0 a <X*X)*I t*n 

TEM11 = TEM8/TEML0 

T EMI 2 = TEM9/TEM10 

TEM13 = PE XP ( X*Tj*OCC$( Y*T ) 

TEM14 = UeXPi X*T) *DSINI Y*T) 

U = <TEM13*TEM11 )-<TEM14*TEM12) 

Y = ( TEM14< t TEMli) ♦ ( TEM13*T£M12) 

RETURN 

END 
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APPENDIX B 



LISTING OF TEST CASES AND EXPERIMENTAL RESULTS 

The Laplace transformation pairs f(s) and F(t) which we 
have used as test cases are listed within this appendix. 

Also included for each function f(s) is the listing of the 
SUBROUTINE VALUE and a typical result of numerical inversion 
at one specific value of t. Case 9 also shows a plot of 
F(t) versus t obtained numerically. 

The results obtained by numerical inversion have been 
compared where applicable, with tabulated values listed in 
(1). In all other cases the value of F(t) has been cal- 
culated on the IBM 360/67. 

The symbols which are used in the presentation of 
results are defined below: 

AL - constant associated with f(s) 

BE - constant associated with f(s) 

T - time for which F(t) is evaluated 

AA - real axis starting position of the integration contour 

x^ - final value of x at termination 

y f - final value of y at termination 

MMMMM - total number of calls made to SUBROUTINE VALUE 

LOSC - number of sign changes between successive addends 

SUM - value of F(t) (equation (16)) 

SUMA - value of G(t) (equation (19)) 

TIME - clock time (sec) required for numerical integration 
ANALYTICAL - analytical value of F(t) 

In all cases we used N = 5 and EPS = l.D-11. 
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Case I 




SU8R0UT I NE V4L'J6 < X , Y , T , AL , U , V , MMmmm j 
THIS IS TABLE ENTRY NUMBER 01 OF L4 PL* 
IMPLICI" P54L *8 ( 4, 3, C-H,0-Z) 

IMPLICIT COMPLEX* 16 (C) 
v(MMMM = + J, 

XPTR = X* T 
XPTI = Y*T 

CXPT = OCM^LX (XPTR, XPTI) 

CNUM = COcXPICXPT) 

CDENOM = OCMPlX(a,Y) 

CRESLT = CMUM/COENO'-l 
U = ORE AL ( CP E$ L T ) 

V = D I MAGI CRESLT) 

RETURN 
E NO 



CE TRANSFORMS BY SPIEGEL 



AL = 1 .253-1 
T = B.D+O 
AA = 0.2 
X F = .4.64 

y? = 2.2 

MMMMM =23 
LOSC = 5 

SUM = 1 . OOOOOD+O 
SUMA = 1.42588D+0 
TIME = O.79794D-I 
ANALYTICAL = 1. OOOOOD+O 
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t 




-l 




SUBROUTINE VALUE ( X, Y , T , A L . U, V , M“ M mm ) 

C T HI 5 IS T ABLE ENTRY N'JM3E= 02 3F LA p LACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL*6 ( A,3,0-H,0-Z) 

IMPLICIT C3MPL EX* 16 (Cl 
MMMMM = MMMMM + l 
C SI = OCMPLXI X, Y) 

XPTR = X *T 
XPTI = Y* T 

C X p T = OCMPLX(X?TR,X p TI) 

CNUM = COEX P( CX 3 T ) 

CDENOM = DCMPLX (X ,Y) 

CDENOM = C S1-CS1 
CRESLT = CN'J U / CDENOM 
U = DAE AL(CRESLT) 

V = DIM AG (CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.2 
Xp - -U . 6^ 

y F • 2.2 

MMMMM = 23 
LOSC = k 

SUM « 8 . OOOOOD+O 
SUMA = 8.04195D+0 
TIME = 0. 532^8D-1 
ANALYTICAL = 8. OOOOOD+O 
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Case 3 




f*- 1 

(*-!)! * 



0 ! = 1 



SUBRO'JT IME VALUE(X,Y , T , 1L , N*OW E* , U, V, MMMMM ) 

C THIS IS TABLE E NTRY NUMBER 03 OF LAPLACE TRANSFORMS 8Y SPIEGEL 

IMPLICIT °EiL^0 ( A,3 ,0-H,0-Z ) 

IMPLICIT C0 M PLEX*16 (C) 

M MM MM = MMMMM + 1 
ZERO = 0.0+0 
CS1 = 0C MD LX (X t Y ) 

XPTR = X*T 
X PT I = Y*~ 

CXPT = DCMPLX (XPTRt XPTI) 

C MUM = COE XP ( C X D T ) 

CDENOM = CSl**NO0WER 
CRESLT = CNUM/COtNOM 
U = ORE AL(CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

ENO 



AL = 1.25D-1 
NPOWER = 3 
T = 8.D+0 
AA = 0.3 
x = -4.54 

F 

y__ = 2.2 
J F 

MMMMM = 23 
LOSC = 3 

SUM = 3 . 20000D+1 
SUMA = 3 • 20008D+1 
TIME = 0.8369^0-1 
ANALYTICAL = 3.20000D+1 
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Case 4 




n > 




r~ i 

l'(n) 



SUBROUTING VALUE Ur Y, AL, POWE 0 , Ur V, NCALU 

C THIS IS TABLE ENTRY NUMBER 04 OF LAPLACE TRANSFORMS eY SPIEGEL 

. IMPLICIT REAL *8 ( A, 3 f C-h , 0-Z ) 

I MPLICIT C0MPLEX*16 ( C) 

DIMENSION NC ALL ( 5 ) 

MMMMM = MMMMM + 1 
CS1 = OCMPLX(X,Y) 

XPTP = X*T 
XFTI = Y*t 

CXPT = DCMPLX(XPTRrXPTI) 

C MU M r CDEXP( CXPT ) 

CALL CPCWER ( CS I r C DENCM r PQW ER r 1 , N CALL ) 

CRESLT = CNUM/COENOM 
I J = DREAL(CPESLT) 

V = 01 MAG(CRESLT) 

RETUPN 

END 



AL = 1.25D-1 
POWER = 1 . 7D+0 
T = 8.D+0 
AA = 0.2 
X p = -4.6U 

y F = 2,2 

MMMMM =23 
LOSC = k 

SUM = U.71815D+0 
SUMA = U.8324D+0 
TIME = 0.85098D-1 
ANALYTICAL = ^.7l8l5D+0 
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Case 5 







SUBROUTINE VALUE (X , Y ,T , AL,'J » V • mmwvm) 

C • THIS IS T*?LE ENTRY NUMBER 05 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL *8 ( A, 3, C-h,0-Z) 

IMPLICIT COMPLEX* 16 ( C) 

MMMMM = + 1 

ZERO = 0.0+0 
CS1 = DCM D LX( X ? Y) 

XPTR = X *T 
XPTI = Y* T 

CXPT = DCMPLXIXPTR, XPTIJ 
CNUM = CDEX PICXPT ) 

CS2 = DC M PLX(ALt ZERO) 

CDENOM = CS1 - CS2 
C RES L T = CNUM/COcNOM 
U = ORE AL(CPESLT) 

V = DIM AG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = -4.99 

y F - 2.3 
MMMMM =23 
LOSC = 5 

SUM = 2.71828D+0 
SUMA = 3 . 59 U 74 D+O 
TIME = 0.89050D+0 
ANALYTICAL = 2. 7182PD+0 



97 



Case 6. 




jn — I e at 

(n- 1)1 



0 ! = 1 



c 



SUBROUT INE VALUE (X,Y,’ f 4Lt NPCWE* ,U, V , MMMMM) 

THIS IS TAB LE ENTRY NUMBER 06 OF LAPLACE TRANSFORMS BY SPIEGEL 
jMPLiri T p E AL *3 ( A, B t 0-F,0-Z ) 

I MPLICIT COMPLEX* 16 (C) 

MMMMM = MM M MM + i 



ZERO = 0.0+0 
C SI = DC MO L X( X » Y) 

CS2 = DCMP LX ( AL t ZERO ) 
XPTR = X*^ 

XPTI = Y*T 

CXPT = DCMPLX(XPT.R f XPTI ) 
CNUM = CDEXP(CXPT) 

C S3 = CS1 - CS2 
CDENOM = CS3**N P 0WER 
C RE SL T = CNUM/ CDENOM 
U = DREAL(CRESLT) 

V = D I M AG ( C RESLT ) 

RETURN 

END 



AL = 1.25D-1 
NPOWER = 3 
T = 8,D+0 
• AA = O.U 

X F * -4.44 

y F - 2.2 

MMMMM = 23 
LOSC = 3 

SUM = 3.69850D+1 
SUMA = 8 , 698650+1 
TIME = 0.7579OD-1 
ANALYTICAL = 8 . 6985 OD+I 
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Case 7. 




**-« e at 
r(n) 



SUBROUTINE VAL'Jc ( X , Y , T . AL , MPOW E R ,U ,V , NC AL L ) 

C THIS IS TABLE ENTRY NUMBER 07 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL *3 ( A, ?, D-l-,0-Z ) 

I MPL I C I " r COMPLEXES (C) 

DIMENSION NC ALL t 5 ) 

MMMMM = MMMMM ♦ 1 

ZERO = C.D+O 
CS1 = OCMCLX(X,Y) 

CS2 = DCMP LX ( AL » ZERO) 

C S3 = CS1 - CS2 
XD-r? = x*7 
X PTI = Y*- 

CXPT = OCMPLXI XPTP., XPTI ) 

CN'IM = CDE X P ( CX °T ) 

CALL CPCWER ( CS3 t CDENCM , ROUE P. , I , NCALL) 

CPESLT = CN UM / CDE NO M 
'J = DREAL(CRESLT) 

V = DI MAG <C RESLT) 

RETJ’N 

END 



AL = 1.25D-1 

POWER = 1.7D+0 

T = 8.D+0 

AA = 0.3 

X =- 4.54 
F 

y F = 2.2 

MMMMM = 23 
LOSC = k 

SUM = 1.28253D+1 
SUMA = 1.30533D+1 
TIME = O. 96928 D-I 
ANALYTICAL = 1.28253D+1 
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Case 8 




3in at 



SU3RTJT IN 6 VALJEI X, Y, T, 41, U, V, MM^m) 

THIS IS T A8LE ENTRY NUM B £R OS OF LAPLACE TRANSFORMS 8Y SPIEGEL 
IMPLICIT REAL*8 ( A ,8 ,0-H ,Q-Z ) 

I ELICIT COM D L EX* 16 ( C) 

MMMMM = MMM«M + I 
Z ERO= 0 . D+O 
CS=OC M PLX ( X , Y ) 

CAL = OCMPLX( 4L ,ZER 3) 

CT=0CMPLX( T.ZERO) 

CDEN = ( <C C *CSH-I CAL*CAU ) 

C6XP=:OEXP( CS*CT) 

C = C EX P/CDEN 
U=OREAL(C) 

V=OI MAGIC) 

RET'JRN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = -U.5U 

y F « 2-2 

MMMMM = 23 
LOSC = ^ 

SUM = 6.73177D+0 
SUMA = 7.03097D+0 
TIME = O. 8629 AD-I 
ANALYTICAL = 6.73177D+0 
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Case 9. 




cos at 



SUBROUTINE VALUE ( X, V, T, &L, U, V, MMMMM) 

C THIS IS t ABLE ENT RY MUM B E° 09 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT p E AL*8 I A,b , D-H, 0-Z) 

I^PL ICIT C0 M P L EX* lo < C ) 

MMMMM = mnvmmm + 1 
ZERO = 0. 0+0 
CS=DCMP LX( X t Y ) 

CAL=0C M PLX ( AL t ZERO) 
c t=o:mplxi T, ZERO) 

CDEN=C3**2+CAL**2 
C=C S/COEN 
C EXP= C oe XP ( CS*CT) 

C =C*C EX P 
U=0RE A L < C ) 

V*0IM AG(C) 

RETURN 
ENO 



“0. 10. 20. 30. 40. SO. SO. 

AL = 1.25D-1 T 

T « 8.D+0 

AA = 0.3 
Xy » -U.99 
Yp - 2.3 
MMMMM = 2U 

Lose = 5 

SUM = 5.^0302D-1 
SUMA * 2 • 00010D+0 
TIME = 0 • 11138D+0 
ANALYTICAL = 5.^0302D-1 
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Case 10 




e' Jt sin at 
a 



SUeRO'JT I ME VAIUE(X,Y ,T f A{_ ,9E, U» V , MMMmM) 

C • THIS IS TABLE ENTRY NUMBER 10 HE LAPLACE TRANSFORMS 8Y SPIEGEL 
IMPLICIT REiL*8 ( A ,8 , 0- H, O-Z ) 

IMPLICIT Complexes <C) 



MMMMM 


— MMMMM + 1 




ZERO = 


0*0 + 0 




CS1 « 


OCMPLX (X, Y) 




CS2 = 


DC MPL X( 3c t Z 


cRO) 


C S3 = 


DCMPLX< AL ,Z 


ERO ) 


X PTR = 


X* T 




XPTI = 


Y*T 




CXPT = 


DCMPLX (X°t 


RtXPTI ) 


C MUM = 


CDE XP { CX PT 


) 


CS4 = 


: SI - CS2 




CDENQM 


= ( CS4*CS4 


)+(C S3* 


C RE $ L T 


= CN ! JM/CO£ 


N3M 


U a OP 


EAL(CRESLT) 




V = 01 


M AG ( C? ESLT ) 





RETURN 

END 



AL = 1.25D-1 
BE = 1.25D-1 
T = 8.D+0 
AA = 0.3 

x p = .4.54 
y F = 2.2 
MMMMM =23 
LCSC = k 

SUM = 1.829B8D+1 
3UMA = 1.88117D+1 
TIME = 0.700l*4D-l 
ANALYTICAL = 1.32988D+1 
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1 



« - 6 

(*'- b)- -r a 2 



c bt cos at 






. SUBROUTINE VALUE! X,Y,T, AL, BE, U,V,w*mmm) 

C THIS IS ^ABLE ENTQY NUMBER 11 OF LA°LACE TRANSFORMS BY SPIEGEL 

IMPLi: IT REAL*3 ( A,6,L>-H,0-Z) 

IMPLICIT C^MPLEX^lo ( C ) 

MNMMM = mmmm.m + i 
ZERO = O.O+O 
CS1 = DCMPLX(X,Y) 

CS2 = OCMPLX! EctZE^O) 

C S3 = OCMPLX(AL,ZEROJ 
X PT R = X*t 
X PTI = Y*T 

CXPT = QCMPLX! XPTR, XPTI ) 

CNUM = (CS 1 -CS 2 ) * CD EXP ( CX p T ) 

C$4 = CS1 - C 52 

CDEMOM = (CS4* CS4) + (CS3^CS3) 

CRES L" = CNIJM/COENOM 
U = OREAL(CRESLT) 

V = OIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
BE = 1.25D-1 
T = 8.D+0 
AA = 0.^ 

X p = -A. 89 
y F = 2.3 
MMMMM = 2U 
LOSC = 5 

SUM = 1.W869D+0 
SUMA = 4.61766D+0 
TIME = 0.7S130D-1 
ANALYTICAL = I.U 6869 D+O 
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Case 12 




sin h at 
a 



SUBROUT INE VALJE(Xf Y» T f AL» Ut V^mmmm ) 

THIS IS TABLE ENTRY MUMBER 12 OP LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT R 5 'L*8 ( A # B, O-HfO-Z) 

I VPLICIT CO M PL EX* 16 ( C) 

MMMMM = MMMMM + 1 

ZEP0=0.D+0 

CS=OCMPLX<X,Y) 

C4L=0CMPLX( AL.2ERQ) 

C’=OC M PLX(T,Z£RD) 

COEN = ( (CS*CS)-< CAL*C4l) 1 

cexp=:dexp<c$*c t> 

C = C EX P/C DEN 
U=QRE4L (C) 

V=DINAG(C) 

RETURN 

END 



AL = 1.25D-1 
T = 9.D+0 
AA = 0.3 

x F « -4.54 
y 7 = 2.2 

MMMMM = 23 
LOSC = 4 

SUM = 9. 40162D+0 
SUMA = 9.64S99D+0 
TIMS = 0.10691D+0 
ANALYTICAL = 9.40162D+0 
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Case 13. 




cosh at 



SUBROUT 
THIS IS 
IMPLICI 
I ^ PLI C I 
MMMMM = 
ZER0 = 0 • 
C$ = DC M P 
CAL=0CM 
CT=DCMP 
COEN = 

c=cs/:o 

CEXP= CD 
C=C*CEX 

u=dreal 

V=DIMAG 
RET J RN 
END 



INE VALJEI X t r'tTt AL, U f V f mmmmm ) 

TABLE ENTRY NUMBER 13 OF LAPLACE TRANSFORMS BY SPIEGEL 
T R E*L*8 ( A, 3 ,D-H,W ) 

T C^mplEX^Io ( Ci 
MMMMM ♦ 1 
D+0 

LX ( X f Y ) 

PLX( AL »Z ER 0) 

LX ( T » Z ER 0 ) 

( (CS^CSl-i CAl*CAL)> 

EN 

EX^ ( CS*CT ) 

P 

(C) 

C C ) 



AL = 1.250-1 
T = 8. D+0 
AA = 0.3 
X p = -4.99 
y F = 2.3 
MMMMM =24 
LOSC = 5 

SUM = 1.54308D+0 
SUMA = 2.52401D+0 
TIME = 0.85176D-1 
ANALYTICAL = 1.5*O08D+0 
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r- v siith at 




SUBPOUT INE V ALi E ( X , Y , T , AL » BE, U, V , MMMMM J 

THIS IS TABLE ENTRY NUMBER 14 CF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT =EAL*8 < A,3 ,D-ri,0-Z) ^ 

I MPLICIT COMPLEX*ib ( C) 

MMMMM = MMMMM + I 
ZERO = O.O+O 
CS1 = DCMPLX ( X f Y ) 

C S2 = OCMPLXt 3E ,ZERO) 

CS3 = DCM° LX ( AL , ZERO J 
XPTR = X*^ 

XPTI = Y*T 

CXPT = CCMPLX(X°TR,XPTI ) 

CNUM = COEXP(CXPTJ 

C $4 = CS1 - C$2 

CDENOM = <CS4*C$4J-(CS3*:$3) 

CRESL" = CNtJM/CDENO A 
U = 3 R E AL (C= ESLT) 

V = DIM AG(CRE$L r ) 

RETURN 

END 



AL = 1.25D-1 
BE = 1.25D-1 
T = 8 . D+0 
AA = 0.5 




MMMMM = 24 
LOSC = 4 



SUM = 2.55562D+1 
SUMA = 2.81918D+1 
TIME = O. 74932 D-I 
ANALYTICAL = 2.555623+1 
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Case 15. 

Z ' 1 



SUBROUT INE V ALUE (X, Y , T , AL »*E » U» V, mmMMM) 

C THIS IS TABLE ENTRY NUMBER 15 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL *8 ( A » 3 f D- ri * O-Z ) 

IMPLICIT CO”PLEX*l6 (CJ 



MMMMM 


a mmmmm + 1 






ZERO = 


0*0+0 






CSI « 


DCMPLX (X » Y ) 






CS2 * 


OC M PL X ( Be f Z 


ERG) 




CS3 a 


DCMPLX(AL, z 


ERO) 




X PTR = 


X* T 






X D T I = 


Y*T 






CXPT = 


DCMPLX (X D T 


Rt X P T I 


) 


CNUM = 


( CSI-CS2 ) * 


CDEXP< 


CXPT J 


CS4 = 


CS1 - CS2 






CDENOM 


= (CS4*C$ + 


J-CCS3 


*CS 3 ) 


C RE SL T 


* cnum/cde 


MOM 




U = DP 


EAL(CPESLT) 






V = 01 


M AG ( CRESL T I 







RETURN 

END 



s — b 



(s — 6)- — a* 



tf bf cosh at 



AL = 1.25D-1 
EE = 1.25D-1 
T = 8.D+0 
AA = 0.5 
X F = -4.79 
Yjr = 2.3 
MMKMM =24 
LOSC = 5 

SUM = 4. I 9 A 53 D+O 
SUMA = 9.63599D+0 
TIME = 0.77740D-1 
ANALYTICAL = 4 . 19453 D+O 
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■1 



1 

(5 - a)(s - 6) 



a * b 



e ht — 
b — a 



£ 



SUBROUTINE VALUE ( X , Y , T , A L , 9E , U f V f v v mmm ) 

THIS IS T ABLE ENTRY NUMBER 16 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT PEAL*3 ( A, 6» C-Ht 3-ZJ 
IMPLICI T COMPLEX* 16 (0) 

MMMMM = MMMMM + 1 
ZERO = 0« D+O 
CS1 = DC MPL X( X » Y) 

CS2 = OCMPLX ( BE » Z ERO ) 

C S3 = DCMPLXULt ZERO) 

XPTR = X*T 
XP~I = Y * T 

CXPT = DCMPLX(XPTR f X°TI) 

CNUM = CDE XP ( CX 3 T ) 

CSA = C$1 - CS3 
C S5 = CS1 - CS2 
CDENOM = CS4*CS5 
CRESL T = CN'JM/CDENOM 
U = OREAL(CPESLT) 

V = DIM AG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
BE = -1.25D-1 
T = 9. D+O 
AA = 0.4 
Xp = -*.89 
y F = 2.3 
MMMMM = 24 
LOSC = 4 

SUM = 9.40161D+0 
SUMA = 1 . 07718D+1 
TIME = 0 . 76648D-1 
ANALYTICAL = 9.40161D+0 
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-1 



(s - a)[s-b) 



a ** b 



be bt — ac at 
6 — a 






SUSROUT INE VALJE(X,Y t 7, AL , 3E, J,V , 

THIS IS TABLE EN T RY .NUMBER 17 OF LAPLACE TRANSFORMS 8Y SPIEGEL 
IMPLICIT REALMS { At 9 f 0-H,0-Z ) 

I VPLICI^ C0MPLEX*16 (C) 

MMMMM = MMMMM + 1 
ZERO = 0.0+0 
CS1 = OCMPLX(X,Y) 

C$2 * DCMPL X( BE ?Z ER 0 ) 

C$3 = DCMPLX { AL, ZERH) 

XPTR = X*t 
XPTI = Y*T 

CXPT = 0CMOLX(X°TR f XPTI J 
CNUM = CSi*CDEXP(CXPT) 

CS4 = CS1 - C S3 
CS5 = CS1 - CC 2 
COENOM = C S4*C S5 
CRESLT = C K| UM /COENOM 
U = OREALfCRESL") 

V = OIMAG(CRESLT) 

RETJRN 

END 



AL = 1.25D-1 
BE = -1.25D-1 
T = 8 .D +0 
AA = 0.3 
X F = .4.99 

y F - 2.3 

MMMMM =24 
LOSC = 5 

SUM = 1.54308D+0 
SUMA = 2.52401D+0 
TIME = O. 78988 D.I 
ANALYTICAL = I. 54308 D+O 
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Case 1 8 




sin at — at cos at 
2a 3 



SLBROUTI NE VALUE (X t Y,T,AL,U,V,MMMmv) 

ThIS IS TABLE ENTRY NUMBER 13 OF LAPLACE TRANSFORMS BY 
IMPLICIT REAL *8 < A, 3» C-Ff O-Z ) 

I MPLIC I T COMPLEX-16 ( C) 

MMMMM = mmmmm + i 
ZERO = 0,0+0 
CS1 = DC MP LX( X * Y) 

XPTR = X*T 
X PTI = Y* T 

CXPT = OCMPLXl XPTR,XPTI ) 

CNUM = COEXP(C XPT ) 

CSA = CSI-CSl 

CS5 = DCMPLXl ALtZERO) 

CS6 = CS5*CS5 
CS7 = CS4 + CS6 
C DEND M = CS7*CS7 
CRESLT = CNUM/CDENOM 
U = ORE AL( CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

ENO 



SPI EGEL 



AL = 1,250-1 
T = 8,0+0 
AA = 0,3 

X y = -4,11 

y F = 2.i 
MMMMM =22 
LOSC = 2 

SUM = 7.709920+1 
SUMA = 7 • 7099^01 
TIME = 0, 51870D-1 
ANALYTICAL 7 • 709920+1 
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Case 19 





( S 2 -r 0 -')2 



t sin at 
2 a 



SUBROUTINE VALUE (X ,Y ,T f ALt U » V ) 

THIS IS TABLE ENTRY NUMBER 19 OF LAPLACE 
IMPLICIT REAL*8 ( A, 9, D-H, O-Z ) 

1 MPLICIT COMPLEX* lo <C) 
mmmmm = MMMMM + 1 
ZERO = 0 • D+0 
CS1 = OC MP L X( X f Y) 

XPTR = X*T 
XPTI = Y*t 

CXPT - GCMPLX (XPTRfXPTI ) 

CS3 = CDEX°(CXPT) 

CNUM = CS1*CS3 

CSA = C S1*C SI 

CS5 = OCMPLXI AL tZERO ) 

CS6 * C S5*C $5 
CS7 = CS4 ♦ CS6 
C CE NO M = C S7*CS 7 
CRESL T = C NUM/ COE NOM 
U = ORE AL(CRESLT) 

V = OIMAGICRESLT J 

RETURN 

END 



TRANSFORMS BY SPIEGEL 



AL - 1.25D-1 
T = 8, D+O 
AA = 0,3 

Xp = -u .54 
yp = 2 • 2 
MMMMM =23 
LOSC = 3 

SUM = 2.69271D+1 
SUMA = 2.69278D+1 
TIME = 0.891S0D-1 
ANALYTICAL 2.69271D+1 
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Case 20 




sin at -f at cos ct 
2 a 



SUBROUTINE VALUE (X * Y »T » A L * 1 J » V 

C THIS IS TABLE E N T R Y NUMBER 20 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL *8 ( A, 3 , C- F , O-Z ) 

I MPLICIT CO MO LE X* 1 6 (C) 

¥ i 

ZERO = 0.0+0 
C$1 = OCMPLX(X,YJ 
C$2 = CSI*0S1 
XPTR =. X* T 
XPTI = Y*T 

CXPT = DCMPLXIXPTR tXPTI J 
CS3 = CDEXO(CXPT) 

CNIM = CS2*CS3 

CS4 = C$1*CS1 

CS5 = DCMP LX ( AL » ZERO ) 

CS6 = C S5*C$5 
CS7 = CS4 + CS 6 
CDENOM = C$7*CS7 
CRESLT = C NUM/C OE NQM 
•J = DREAL (CPESLT) 

V = DI MAG (CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = . 4.54 
yp = 2.2 
MMMMM = 23 
LOSC = 4 

SUM = 5.52709D+0 
SUMA = 5.89548D+0 
TIME = 0.60634D-1 
ANALYTICAL 5. 52709D+0 
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Case 21 




S 3 

(s* + 



cos at — ^ at sin at 



SUBROUTINE VALUE ( X, Y, T, 1L , U, V, MMMMM) 

C THIS 1$ T A9L E ENTRY NUMBER 21 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL*8 < A,3 ,D-H,0-Z) 

IMPLICIT CO M P L EX* 16 ( C ) 



MMMMM 


— M M M M M + i 


ZERO 


= 0* D+ 0 


CS1 = 


DC MP LX ( X , Y ) 


CS1 s 


DC M0 LX ( X f Y ) 


C $2 = 


C $l*CSl*C3i 


X PTR 


= X*T 


XPTI 


= Y*T 


CXPT 


= DCMPLX( XPTR, XP 


C$3 = 


CDEX P(CXPT) 


CNUM 


= C$2* C S3 


CS4 = 


CS1*CS 1 


CS5 = 


DC M P LX ( AL » ZERO) 


C 56 = 


CS5*C$5 


CS7 = 


C$4 + CS6 



CDEMOM = C*7*C$7 
CRESLT = CNUM/COENOM 
U = DR E£L( CRESLT) 

V = DIMA61CRESLTJ 

RETURN 

ENO 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = -4.99 

yp = 2.3 

MMMMM = 24 
LOSC = 5 

SUM = 1.19567D-1 
SUMA = 1.89481D+0 
TIME = 0.10865D+0 
ANALYTICAL 1.195670-1 
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Case 22 




t cos at 



SUB*1U T IME V^LUE (X ,Y f T, AL,U,V,MMMM m) 

this IS TABLE ENTRY NUMBER 22 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT RE 4L *8 ( A,B,0-H,0-Z) 

IMPLICIT C0«PLEX*l6 ( C J 
MMMMM = MMMMM + 1 
ZERO = O.D+O 
CS1 = DC MD LX (X , Y ) 

C'2 = OCMPLX(AL.ZERO) 

XPTR = X*T 
X PT I = Y *T 

CXPT = DCMDtX(XPTR,XPTI ) 

CMUM = (CCS 1-C S 1 J - (CS2*CS2> ) *C D E XP( C X PT ) 

CR3 = (<CS1*CS1I + (CS2*CS2)J 
CDENOM = C$3*CS3 
CRESL T = CNUM/COENOM 
U = DREAL(CRESLT) 

V = DIMAG(CRESLT) 

RETJRM 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
Xp = -n .54 
yp = 2.2 
MMMMM =23 
LOSC = 4 

SUM = 4 . 3224210+0 
SUF1A = 4 . 76 OOOD+O 
TIMS = 0.6721CD-1 
ANALYTICAL 4.32242D+0 
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-1 



1 

(« 2 -a 2)2 



at cosh at — sinh a£ 
2a* 



cF 



SUBROUTINE VALUE (X , V , T OL, J , V , ) 

THIS IS TABLE ENTRY NUMBER 23 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT 0 EM *8 ( A,B,0-H t 0-2) 

I VPLICIT C3 M PLE:< J ‘I6 ( C ) 

MMMMM = MMMMM + I 
ZEPO = 0,0+0 
C$1 = 0CMPLX(X,Y) 

XP TP = X*T 
X FT I = Y *7 

CXPT = DCMPLX {XPTR t XPTH 
CNUM = COE XP( CXPT) 

C$4 = CSl*C$i 
C $5 = DCMP LX ( AL ,Z£RQ) 

CS6 = CS5*CS5 
CS7 = CS4 - CS6 
CDENOM = C S7*CS7 
CRESL T = CNUM/COENJM 
U = 0 RE AL (C RE$L T ) 

V = 01 MAG ( CRESLT) 

RET J RN 
ENO 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4.44 

y F = 2.2 
mmmmm = 23 
losc = 2 

SUM = 9.^17712+1 
SUMA = 9.^17712+1 
TIME = 0.61984D-1 
ANALYTICAL 9.^1771D+1 
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Case 24 




t sinh at 
2a 



C 



S LB 5 OUT I ME VALUE l X , Y f T , A L. U » V 

THIS IS T A 5 L E E'\T R Y NUMBER 24- 

IMPLICIT REAL*8 l A , S , C-H , >Z ) 

IMPLICIT CHMPLEX*16 1C) 

mmmmm = MMMMM + 1 

ZERO = O.O+O 

CSI = DCM^L X( X » Y ) 

XPTR = X* T 
XPTI = Y*t 

CX°T = OCM°LX(XPTR»X D TI ) 

C S3 = CDEX P( CXP T ) 

CNUM = CSI*CS3 

CS4 = CSI*CS1 

CS5 = DC M PLX(AL,ZERO) 

CS6 = C S 5*C S 5 
CS7 = CS4 - C56 
CDENDM = C$7*CS7 
C RE S L T = CNUM /C Dc NOM 
U = ORE AL (CRESLT ) 

V = 01 MAG(CPESLT) 

RE T JRM 

END 



9 MMMMM) 

OF LAPLACE TRANSFORMS BY SPI EGEL 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4.44 

y F = 2.2 
MMKMM =23 
LOSC = 3 

SUM = 3 • 76064D+1 
SUMA = 3.76151D+1 
TIME = 0.66924D-1 
ANALYTICAL 3 . 76064D+1 
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Case 25 




sinh at + at cosh at 
2a 



SUBROUTINE VALUE (X, Y,T, AL,U,V,MMMMM) 

THIS IS TABLE ENTRY NUMBER 25 OF LAPLACE TRANSFQP M S BY SPIEGEL 
IMPLICIT REALMS ( A, 8» D-H, 3-Z ) 

IVPL1CIT C0MPLEX*16 (C) 
mmmmm = M^MMM + 1 

ZERO = O* D + 0 
CS1 = OCYPLX (X i Y ) 

CS2 = CSI-CSl 
XPTR = X* T 
XP T I = y*T 

CXPT = DCMPLX(XPTR,XPTn 
CS3 = CDEX °( C XP T ) 

C NJM = CS2*CS3 

CS4 = CSl*CSl 

CSS = DCMP LX ( AL t Z ERO ) 

CS6 = C $5*CS5 
C $7 = CS4 CS6 
CDENO* = C S7-CS7 
CRESL T = CNIJM/COENOM 
U = 3PEAL(CRESL T ) 

V = 01 M AG( C PESLT ) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = O.k 




MMMMM = 2 1 * 
LOSC = U 



SUM = 1.08731D+1 
SUMA = 1.20870D+1 
TIME = 0.73710D-1 
ANALYTICAL = 1.08731D+1 
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Case 26 




cosh ai + smhat 



SUBROUTINE VALUE ( X , Y , T , A l,'J , V , mmmmm) 
THIS IS T A RLE ENTRY NU*3ER 26 OF LAPLACE 
IMPLICIT REAL-8 ( A , 8 , C-h , O-Z J 
IMPLICIT COMPLEX* 16 (C) 

MMMMM = MMMMM + i 
ZERO = O.O+O 
C SI = DC MP L X( X, Y) 

CS2 = CSI*CS1*CS1 
XPTR ^ x* T 
XPTI = Y*T 

CXPT = OCMPLX (X^T R, XPTI ) 

C S3 = CDEXP(CXPT) 

C MUM = CS 2*C S 3 

CS4 = CS1-C51 

CS5 = DC M PLX( AL, ZERO) 

CS6 = CS5*CS5 
CS7 = CS4 - C36 
C0EM3M = C S7-C S 7 
CRESL T = C\«;jM/COENOM 
U = ORE AL(CRESLT) 

V = OIMAG(CRESLT) 
retjrn 

END 



TRANSFORMS BY SPIEGEL 



AL = 1.25D-1 
T = 8 . D+0 
AA = 0.4 
X F = -4.89 
y F = 2.3 
MMMMM = 24 
LCSC = 5 

SUM = 2 . 13068D+0 
SUMA = 4.30908D+0 
TIME = 0.19490D+0 
ANALYTICAL 2. 13068D+0 
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Case 27 




t cosh a t 



SUBROUTINE VALUE ( X f Y,T,AL»U*V*MVMMM) 

C THIS IS TABLE ENTRY NUMBER 27 OF LA°LAC E TRANSFORMS BY SPIEGEL 

IMPLICIT REAL*8 { A t 3 .D-h,0-Z) 

IMPLICIT COMPL EX* 16 (C ) 

MMMMM = MMMMM + I 
ZERO = O.D+O 
CS1 = DCM D L X( X, Y) 

CS2 = OCMPLX (AL, ZERO) 

XPTR = X*T 
X D T I = Y*T 

CXPT = CCMPLX(XP^R,XPTI ) 

CNUM = <(CSl*CSl) + (C$2*CS2))*C0EXP(CXPT) 

CS3 = <(C$1*C51) - ( CS 2*CS2) ) 

C DENQM = CS3*CS3 
CRESLT = CNUM/COENOM 
U = OR E AL ( r R ESLT ) 

V = 0 I M AG ( C RESLT ) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
- AA = 0.4 
X F = -4.89 
y-p - 2.3 
MMMMM =24 
LOSC = 4 

SUM = 1.23446D+1 
SUMA = 1.34023D+1 
TIME = 0.65078D-1 
ANALYTICAL 1.23446D+1 
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Case 28 




(.1 — a 2 t~) sin rtf — 3rtfcosa< 

8<i i 



SUBROUTINE VALUE (X,Y,T»ALfU»V,MMMMV) 

THIS IS ^ A 8L E ENTRY NUMBER 29 OF LAPLACE TRANSFORMS BY SPIEGEL 
I MPLIC I T REALMS ( A, 8, C-htO-Z) 

IMPLICIT C0MP L ex^l6 (C) 



MMMMM 


s MMMMM + l 


ZERO 


= 0 • D+0 


CS1 = 


DCMPLXC X t Y) 


X PTR 


= X*T 


X PTI 


= y* t 


CXPT 


= DCMPLXt XPTR, X° TI ) 


CMJM 


= CDEXP(CXPT) 


CS4 = 


CS1*CS1 


C S5 = 


DCMPLXC AL tZERO ) 


CS6 = 


CS5*CS5 


CS7 = 


CS4 * CS6 



COENOm = CS7*CS7*CS7 
CRESLT = CNUM/COENOM 
U = 0 REAL (CRESLT) 

V = DIM AG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
Xjr = -3.6 

yp = 2.o 

MMMMM = 21 
LCSC =■ 2 

SUM = 2 . 54096 D+ 2 
SUMA = 2.89148D+2 
TIME = 0.50414D-1 
ANALYTICAL = 2 . 54096D+2 
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Case 29 




t sin at — n t : cos at 
8a J 



SUBROUTINE VALUE (X , Y , T , A L , U , V f vmm VV ) 

C THIS IS T^BLE ENTRY NUMBER 29 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT RE£L*3 ( a,B, D-I-.O-Z ) 

IMPLICIT C0MPLEX*16 ( C) 

MMMMM = MV»*J|NAM + l 

ZERO = 0-0+0 
CS1 = 0CMPLX(X,Y) 

XPTR = X* T 
XPTI = Y*t 

CXPT = DCMDLX(XPTR,XPTI ) 

CS3 = CDEX D ( CX ) 

CMJM = CS1*CS3 

CS4 = CS1*0S1 

CS5 = 0CMPLX(4L, ZERO) 

CS6 = C$5*CS5 
CS7 = CS4 + CS6 
COE MOM = CS 7*CS7* CS7 
CRESLT = CNUM/COENOM 
U = 0° E AL ( C 3 E SL T ) 

V = OIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = -4.11 
y F = 2.1 
MMMMM = 22 
LCSC = 3 

SUM = 1.54198D+2 
SUMA = 1.66759D+2 
TIME = 0.56472D-1 
ANALYTICAL = 1.5^1990+2 
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Case 30. 



'1 




8- 

(*- + a -) 3 



(I + o-V-) sin at — at cos at 



SUBROUTINE VALUE { X , Y, T ,AL , U, V , 

THIS IS T A BL E ENTRY NUMBER 30 OF LAPLACE TP AN SFQP.MS BY SPIEGEL 
IMPLICIT REAL*8 ( 4,8,0-H f 0-Z) 

IMPLICIT COMPLEX* 16 < C ) 
mmmmm = MMM.MM + 1 
ZERO = 0.0+0 
CSI = DC M PL X ( Xf Y ) 

CS2 = CS1-CS1 
XPTR = X*T 

XOTj = y*T 

CXPT = DCMPLX (XPTR,XPTI) 

C S3 = COEXP(CXPT) 

CMJM = CS2-CS3 

CS4 » CSI*CSI 

CS5 = 0C MO LX ( AL y Z ERO ) 

CS6 = CS 5*CS 5 
CS7 = CS4 + CS6 
CDENOM = C $ 7*C S 7*C S7 
CRESL T = CNUM/COENOM 
U = OREALtCRE SLT) 

V = DIM AG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = -4.11 

y F = 2.1 

MMMMM = 22 
LOSC = 2 

SUM = 7.312893+1 
SUMA = 7.31289D+1 
TIME = 0.57122D-1 
ANALYTICAL = 7-31239D+1 
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Case 31 • 




3 / 3 m at f of- cos at 
Ha 



SUBROUTINE VAL J E ( X , Y , T , A L, U, V, mmmmM) 

C THIS IS * r AR LE E NT RY NUMBER 3i OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT ^EAL-3 ( A,B ? D-H,3-Z) 

IMPLICIT COMPLEX* 16 ( C) 

MMMMM = MMMMM + 1 

ZERO = 0.0+0 
CS1 * DC MO LX ( Xf Y ) 

CS2 = CSl*CSi*C$l 
X°TR = X*T 
XP T I = Y*^ 

CXPT = OCMPLXI XPTR, XPTI } 

C S3 = CDEXP( r X°T) 

CNUM - CS2*C$3 

C$4 = CSI*CS1 

CSS = DCM3 LXlALt ZERO) 

C $6 = C$5*CS5 
CS7 = CS4 + C So 
COEMOM = C^ 7*C S 7* CS 7 
CRESL^ = CNLJM/CucNOM 
U = OREAL(CPFSL^) 

V = OIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
x F = ,4.54 
y F = 2.2 
MMMMM =23 
LOSC = 3 

SUM = 2.45177D-I-1 
SUMA = 2.451640+1 
TIME = 0.605540-1 
ANALYTICAL = 2.451770+I 
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Case 32. 




(3 — q J /-) sin at 4- 5at cosat 
Sa 



SLBR3UTINE VA LUE ( X , Y ,7, ALt U t V , 

C THIS IS TABLE EnITRY NUMBER 32 3F LAPLACE TRANSFORMS BY SPIEGEL 

I MPLIC I T REALMS ( A, 3t C-htO-2) 

IMPLICIT CO^PLEX*lo (C) 

MVMMM = + I 



ZERO 


= 0.0+0 


CSX = 


DC MPL X ( Xf Y) 


CS2 = 


CS 1*CS1*CS1*CS1 


X PTR 


= X* T 


XPTI 


= Y*T 


CXPT 


= DCMOtX (XP7R,XPT I) 


C S3 = 


CDEXP(CXPT) 


C NUM 


= CS 2*0S3 


C54 = 


CS1*C$1 


C$5 = 


OCMP LX( A L t ZE R3 J 


CS6 = 


CS5*C$5 


CS7 = 


CS4 + • CS6 



Uu ' » J 'M * ^ O I w J I ' U V ( 

CRESL T = CN'JM/COtNOM 
U = D RE AL(0RESL' r ) 

V = DIMAG(CRESLT) 

RE T JRN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = -u .54 

y F = 2.2 
MMMMM = 23 
LOSC = 4 

SUM = 4.38445D+0 
SUMA = 4.83048D+0 
TIIAE = 6.7548D-2 
ANALYTICAL = 4.384450+0 
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Case 33 




(8 — a 2 t-> cos af — lat sin at 
8 



SUBROUTINE VALUE IX ,Y ,7 , A Lt U , V , ) 

C THIS 1$ T&9LE ENTRY NUMBER 33 OF LAPLACE TRANSFORMS BY SPIEGEL 

implicit realms i a, 9, d-h, o-z ) 

I VPLICIT COMPLEX* 16 (C) 



MMMMM 


= MMMMM + 


i 


ZERO = 


0 • D+0 




CSI = 


DC M P L X( X, 


Y) 


CS2 = 


: si*csi*C 


Sl*CSi*CS 


XPTR = 


X*T 




X PTI = 


Y*T 




CXPT = 


DCMOLX(XPTR*XPTI ) 


CS3 = 


CDEX 0 (CXP 


• ) 


CNUM = 


C S2*C S3 




CS4 = 


CS1*CS1 




CS5 = 


DCMPLX ( AL 


, ZERO) 


CS6 = 


C $5*C S5 




CS7 = 


CS4 ♦ CS6 




CDENOM 


= CS7*CS 


7*C$7 


CRESLT 


= :num/coenom 



V = DI MAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 

Xp = -4.99, 

y F = 2.3 

MMMMM =24 

Lose = 5 

SUM = -2.63523D-1 
SUMA = 1.81943D+0 
TIME = 6.74414D-2 
ANALYTICAL = -2.63523D-1 
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t 2 sin at 
2 a 




3a 2 - a 2 
(s 2 + a-) 3 



SUBROUTINE VALUE (X ,Y ,T » AL,U , V , mmmmm ) 

THIS IS TABLE ENTRY NUMBER 34 OF LAPLACE 
IMPLICIT REALMS { A t 8, D-h, l-l) 

I MPLICI T C Q M0 L£X* 16 ( C i 
MMMMM = MMMMM ♦ 1 

ZERO = O.D +0 
CSl = DC MPL X ( Xt Y) 

CS2 = OC^o L X(AL ,ZERO ) 

X PTR = X* T 
X° TI = Y*T 

CXPT = OC M PLX(XPTR,XPTn 

CNUM = ( (3.D+0*(CSl*C$m-(C$2*rS2) )*COEXP(CXPT ) 
C S3 = <CSI*C$1) <CS2*CS2) 

CDENOM = CS 3=*C$ 3* CS 3 
CRESLT = CN1JM/CDENOM 
U = D® EAL(CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

END 



TRANSFORMS BY SPIEGEL 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 

x F = -k.hh 
y F = 2.2 

MMMMM = 23 
LOSC = 2 

SUM = 2.15417D+2 
SUMA = 2.15U17D+2 
TIME = 0.67834D-1 
ANALYTICAL = 2.151^7D+2 
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Case 35 




8 3 — 3rt 2 s 
(• 2 + 



£f 2 cos at 



SUBROUT INE V&LUE ( X , Y , T.A L,U,V , MMMMM) 

C T HIS IS TABLE entry number 35 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICI T RE&L*8 ( A ,8 ,D-H f G-Z ) 

IMPLICIT C0MPL£X*16 1C) 

MMMMM = MMMMM + 1 
ZERO = 0.3+0 
CS1 = DCMDLX(X,YJ 
C$2 = DCMPLXI AL, ZERO ) 

XPTR = X*T 
X PT I = Y*T 

CXPT = DC M P LX ( X P T R » X °T I ) 

CNUM = ( (C 31*CS 1*CS1) -< 2. 3+0*CSl*C$2*C$2 ) )*CDEXP(CXPT) 

CS3 = (CS1*CS1) + (CS2*CS2) 

CDENOM = CS 3=CS3* CS3 
CRESLT = CNUM/CDENO'-I 
'J = ORE Al (CFESL T ) 

V = 01 MAG(CSESLT) 

RETURN 

ENO 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 

x F = -k.5k 
y F = 2.2 

MMMMM =23 
LOSC = 3 

SUM = 1.72897D+1 
SUMA = 1.72904D+1 
TIME = 8.AA7AD-2 
ANALYTICAL = 1.72S97D+1 
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Case 36 




s 4 — Ga-s- - 1 - a* 
(a* + a -) 4 



£t 3 cos at 



SUBROUTINE VALUE (X,Y,T,4L,U,V,MMM^m) 

C THIS IS TABLE ENT R Y NUMBER 36 LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REALMS ( A, Bf C-H t O-Z) 

IMPLICI T COMPLEX- 16 (Ci 
MVMMM = + 1 

ZERO = O.O+O 
CS1 = DC M D L X< X f Y) 

CS2 = DCMPLX ( ALt ZERO) 

C S3 = CS1*C$1*C$1*C^ 1 
CS4 = 6.D + 0*C$2*C $2*C$1*C$1 
C$5 = CS2*C$2*CS2*C$2 
X PTR « X*T 
XPTI = Y*T 

CXPT = OCMPLX (XP t R»XP t I ) 

CNUM = (CS3-CS4+C S5)-CDEXP(CXP T ) 

CS6 = ((CSI-CS1) + ( C 3 2*CS 2 ) ) 

CDENOM = C56*C5o*C$6*C56 
CRESLT » CNUM /CDENOM 
U * DPEAL(CPESLT) 

V * DI MAG < C RESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA =0.4 
X F = -4.44 

y F - 2.2 

MMMMM = 23 
LOSC = 2 

SUM = 4.61058D+1 
SUMA = 4.61058D+1 
TIME = 0.73424D-1 
ANALYTICAL = 4.61058D+1 
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Case 37- 




< 3 sin at 
24a 



SUBROUTINE VALUE ( X , Y , T f A L , U , V , MM mmm) 

C THIS IS TABLE E47RY NUMBER 37 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL* 3 ( A ,3 * D-H » C-Z ) 

IMPLICIT r -0 *°LE X* 16 (C) 

MMMmm = MMMMM ♦ i 

ZERO = C.0+0 
CSI = DCM?LX(X,Y) 

CS2 = DCMPLXCAL, ZERO) 

XPTR = X*T 

XPTI = Y*T 

CXPT = OC M PLX(XPTq f xP T I) 

CNUM = < (C Sl*CSi*C$I) - (CS2*CS2*C$1 ) ) *CDE XP ( CX PT ) 

CS3 = UCS1*C51) + ( C$ 2 *CS2 ) ) 

C OE NO M = CS3*CS3*C$3*C$3 
CRESLT = CNJN/CJENOM 
U = 0 RE AL ( C D £$L T ) 

V « 01 M AG( CR.ESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X p = -4.01 
y F = 2.1 
MMMMM =22 
LOSC = 3 

SUM = 1.43611D+2 
SUMA =1.44005D+2 
TIME = 0.68692D-1 
ANALYTICAL = 1.43611D+2 
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1 

(#2-a2)3 



(3 + g2f 2 ) s inh at — 3a< coshat 
8a* 




S U8RGU T INE VALUE (X*YfT» AL * U» V» mmmmm ) 

C THIS IS TABLt t v|T RY NUMBER 38 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT PEAL-3 ( A, R , D-H ,0-Z ) 

I VPLICIT COM PL EX* lo ( C) 



MMMMM 


- MMMMM + 1 


ZE°0 = 


0.0 + 0 


CS1 = 


DCMPLXIX, Yi 


XP TR = 


X^T 


XPTI = 


Y*T 


CXPT = 


CCMPLX (X?TR,XPT I ) 


C NUM = 


CDE XP(CXPT) 


CS4 = 


CSl*CS 1 


CS5 « 


OCMPLX < AL, ZERO) 


CS6 = 


C S5*CS5 


CS7 = 


CS4 - CS6 


CDENOM 


= CS7-CS7-CS7 


C&ESLT 


= CNUM/COENOM 


U = OR 


EAL (CRESLT j 


V = 01 


MAGCCRESLT) 


RETJRN 




END 





AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X p = -4.01 
y F = 2.1 
MMMMM = 22 
LOSC = 2 

SUM = 2.93122D+2 
SUMA = 3.6?4?lD+2 
TIME = 0.48334D-1 
ANALYTICAL = 2.93122D+2 
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Case 39 




ct 2 cosh at — t sinh at 
8a 3 



SUBROUTINE VALUE (Xf Y,T, ALtU, V ? vmmmm) 

C THIS IS "’ABLE ENTRY NUMBER 39 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT RE^L-3 ( A,BfD-H,a-Z) 

1 PPL I C I T C^hdlex^ 16 (C) 

MRMMM = MMMMM + 1 
ZERO = C.D+O 
CS1 = OCM°L X ( X t Y ) 

X PTR = X*T 
XPTI = Y*T 

CXPT = DCM D LX (X PT R * X PT I ) 

C S3 = CDEXP(CXPT) 

CNUM = CS1-CS3 

CS4 * CS1*C51 

CS5 = OCMPLX(AL» ZERO! 

CS6 * CS5*C^ 5 
CS7 = CSA - CS6 
CDEN3M = CS7*CS7*CS7 
CRESL" = CNUM/COENOM 
U = ORE AL(CRESLT) 

V = OIMAG(CRESLT) 

RE T URN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4.01 

y F = 2.1 

MMMMM =22 
LCSG = 3 

SUM = 1.88354D+2 
SUMA = 1.891953+2 
TIME = 0.48750D-1 
ANALYTICAL = 1.38354D+2 
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■1 



(g2 - 



at cosh at -+ ( a 2 t - — 1) sinh at 
8a 3 



£ 



SLBROUTINE VALUE ( X , Y , T , A L , U , V , VMM v m ) 

THIS IS TABLE ENTRY NUMBER ^0 OF LAPLACE TRANSFORMS 6Y SPIEGEL 
IMPLICIT REAL-6 (AfBfC-HtO-Z) 

implici t :omplex* i6 <o 

MHMMM = m.^mm + 1 
ZERO = 0.0+0 
c SI - dcmplxi X, Y) 

CS2 » C$i*CSl 
XPTR = 

XPTI = Y*T 

CXPT = DCM°LX (X PT R t X°T I ) 

C S3 = CDEXF(CXPT) 

CNU'I = CS2*CS3 

CS4 = CS1*CSI 

CS5 = DCMPLX(AL,Z5R3) 

CS6 = C$5*CS5 
C$7 = CS4 - CS6 
COENOM = C ^7*CS7*CS7 
CP E$L T = CN'JM/COENIM 
U = OREAL(CRESLT) 

V = DIMAG(CPESLT) 

RETURN 

ENO 



AL = 1.25D-1 
T = 8 . D+0 
AA = 0.4 
X F = -U.44 
y F = 2*2 
MMMMI4 =23 
LOSC = 2 

SUM = 9.87572D+-1 
SUMA = 9 • 37572D+1 

time = 0 . 5252 D -1 

ANALYTICAL = 9.37572D+I 
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Case ^1 




3tsinhat 4- at 2 coshat 
8a 



C THIS IS ^ABLE £N T R Y NUMBER 41 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL-8 < A , 3 , D-H , O-Z) 

implicit cqmplex*16 (cj 



MMMMM 


= MMMMM 4 


I 


ZERO = 


C.D+O 




CS1 = 


OCMDLX (X, 


Y) 


CS2 = 


CS1*C.S1*C 


SI 


X D TR = 


X*T 




XPTI = 


Y*T 




CXPT = 


OCMOLXIX 


P T R » X P 


CS3 = 


CDEXP ( C XP 


T) 


CNJM = 


C^2* CS 3 




CS4 = 


C S1*C SI 




CS5 = 


OCMDLX ( AL 


,ZERO ) 


CS6 = 


C $5 *C 55 




C$7 = 


CS4 - CS6 





CDENOM = C S 7*CS 7* CS 7 
CRESL T = CNUM/COENOM 
U = DREAL(CRESlT) 

V = D I M AG ( C RESLT ) 

RETURN 

ENO 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
Xp = -4.44 
= 2.2 

MMMMM = 23 
Lose = 3 

GUM = 4.05495D+1 
SUMA = 4.05582D+1 
TIMS = 0.53508D-1 
ANALYTICAL = 4 . 05495 D+I 
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Case 42 




(3 + a 2 f 2 ) sinh at -f 5at cosh at 
8a 



SUBROUTINE VALUE < X ,Y ,7 ,AL,’J , V , MMVmm) 

C THIS 13 "iBLE ENTRY NUMBER 42 3F LAPLACE TRAM SFHRM S BY SPIEGEL 

IMPLICIT REALMS ( A f 6 f 0-n , 0-Z ) 

IMPLIC I T COMPLEX- 16 < CJ 
HMMMM = mmmmm + x 
ZERO = 0*0+0 
C SI ■ OCMPLXC Xt Y) 

CS2 = CSl*CSl*C$i*CS 1 

XPTR = 

XPTI = y*t 

CXPT = DC MO LX (X DT R# XPT I ) 

C S3 = CDEXP(CXPT) 

CNUM = CS2-CS3 

CS4 = CS1*CS1 

CS5 = OCMPLXt AL f ZERO) 

CS6 = CS5*CS5 
CS7 = CS4 - C 36 
CDEN3M = C S7-C S 7*C S 7 
CRESLT = CNUM/COENQM 
U = DREAUCRESLTJ , 

V = OIMAG(CRESLT) 

RE T 'JRN 

END 



AL = 1.25D-1 
T = 8 .D +0 
AA = 0.4 
Xp = -A . 89 

y-p ~ 2.3 

MMMMM = 24 
LOSC = 4 

SUM = 1.24162D+1 
SUMA = 1.34824D+1 
TIME = 0.582U0D-1 
ANALYTICAL = 1.24162D+1 
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Case k3 




(8 + a-£ 2 ) cosh at + 7afsinha* 
8 



SUBROUTINE VALUE (X ,Y ,T , AL,U ,V , MMMMM ) 

C THIS IS TABLE ENTRY NUMBER 43 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT PEAL *8 ( A, 3 , C-h,0-Z J 

I MPLICIT Cn*OLEX*16 ( C) 



MMMMM 




1 




ZERO 


= 0 .0+0 






CS1 = 


DCMPLXtX, 


Y) 




CS2 = 


C$l*C$i*C 


Si*C 


Si *cs 


X FT R 


= X*T 






XPTI 


= Y* T 






CXPT 


= CCM PL X ( X 


D TR , 


XPTI ) 


CS3 = 


C DEX D ( CX P 


7) 




C NUM 


= CS2*CS3 






CS4 = 


C S I * r S 1 






C$5 = 


DCMP LX ( AL 


♦ ZERO ) 


CS6 = 


C S5*C S5 






CS7 = 


CS4 - CS 6 







CDEN3M = CS7*CS7*CS7 
CRESLT = C NU M /C Ot NOM 
U = ORE AL (CRESLT ) 

V = 01 MAC (CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4.89 
y F - 2.3 
MMMMM =24 
LCSC = 5 



SUM = 2.76427D+0 
SUMA = 4.73261D+0 
TIME = 0.62920D-1 
ANALYTICAL = 2.76427D+0 
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Case kk . 




t 2 sinh at 
2 a 



SUBROUTINE VALUE IX ,Y ,Tf4L,UfVf 

THIS IS TABLE ENTRY NUMBER 44 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICI T ?E*L*3 ( A, 6 f D-H, >Z ) 

IMPLICIT COMPLEX^ 16 ( C J 
MMMMM = MMMMM + 1 
ZERO = O.D+O 
CS1 = OCMPLX(X,Y) 

CS2 = DCMPLXUL,7.ERO) 

XPTR = X*’ 

XPTI = Y 

C. XPT = OCMPLX(X°TR,X°TI ) 

CNJM = t (3 .0+0* (CS 1-CS 1) ) +(CS2*CS2 ) ) *COE XP { CXPT ) 

C S3 = (CS1*CS1) - ( CS2 *C S2 ) 

CDEMOM = CS3*'S3*CS 3 
C RES L T = CNU^/COENT* 

U = ORE At(CRESlT) 

V = 0 I M AG( CRESLT ) 

RE T URN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = O.U 
Xp = -4.4U 

y F = 2.2 

MMMMM =23 
LOSC = 2 

SUM = 3 • 00852D+2 
SUMA = 3.00852D+2 
TIMS = 0.55042D-1 
ANALYTICAL = 3.OO052D+2 
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Case 45. 




cosh at 



SL8R0UT I NE VALUE (X ,Y ,T , A Lt U » V , VMM ) 

C THIS IS TABLE ENTRY NUMBER 45 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT 3 EAL*3 < A, 8, C- h, .0- Z ) 

I MPLI Cl T COMPLEX- Io ( CJ 
MMMNM = MMMMM + 1 
ZERO = 0.D+0 
CS1 = DC L X( X» Y) 

CS2 = DC MP L X( Al t Z ER 3 ) 

XPTR = X*T 
X PTI =• Y*T 

CX°T = DCMDLX(XPTR. XPTI) 

C NJ M “ ( (CSl*C$l*C*i)+(3 .0+0 S1*CS2*CS2I ) *CDE XP { C X°T ) 

C S3 = (CSI*C$1) - (CS2*CS2) 

CDENOM = CS3*CS3*:S3 
CRES L T = C N'JM/ CDENOM 
U = DREAL(CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

END 



AL =■ 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4.44 
y F = 2.2 
MMMMM = 23 
LOSC = 3 

SUM = 4.93786D+1 
SUMA = 4.93892D+1 
TIME = 0. 57772D-1 
ANALYTICAL =• 4.93786D+1 
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Case 46 





+ n‘ 
«-)* 



£ f 3 cosh at 



SU8R3UTINE VALUE ( X ♦ Y , T , AL,'J , V , MMMmm } 

C THIS IS TABLE EN“RY NUMBER 4t> OF LAPLACE TRANSFORMS BY SPIEGEL 

HPL IC IT Q EAL*3 ( AfB»0-H»0-Z) 

I MPL I C I T CO Mp LEX* 1 6 (Ci 
MMMMH = + 1 

ZERO = 0*0+0 
CSI = 0CMPLX(X,Y) 

CS2 = DCMPLXI AL,ZER0) 

CS3 = CS1*CSI*CS1 *CS 1 
CS4 = 6# 0+0 *CS2 *C S2 *C SI *CS1 
CSS = C S2*CS2*CS2*C$ 2 
X pr R = X*T 
XPTI = Y*T 

CXPT = OCMPLXIXPTR, XPTI ) 

CNJM = (C$3 +CS4+CS5 ) *COEXP( CX D T ) 

CS6 = (<C$l*C$l) - (C$2*C$2)) 

CDENOM = CS6*CS6*C$6*CS6 
C RESl T = CN1JM/C06NQM 
U = DPEAl(CRESLT) 

V = 0 I** 4G( CPESLT ) 

RETURN 

END 



AL = 1.25D-X 
T = 8.D+0 
AA = O.k 
X p = 

y f = 2.2 

MMMMM = 23 
LOSC = 2 

SUM = 1.31676D+2 
SUMA = 1.31676D+2 
TIME = O. 65676 D-I 
ANALYTICAL = I. 31676 D +2 
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Case k7 




a- 1 4- a'** 
^ rrr ^¥ 




t 1 sinh at 
21 a 



SLBRO UTINE VALUE (X , Y ,T # AL»*J • V . MMMMM) 

THIS IS ~ A B LE E 'nIT RY NUMBER 47 JF LAPLACE TRANSFORMS 8Y SPIEGEL 
IMPLICIT REALMS < A, 3 t C-Ht H-Z) 

I M°L I C I T COMPLEX- 16 ( C) 

M m MMM = + 1 

ZERO = O.D+O 
C SI = DCMPLX( x, yj 
CS2 = DC MP LX ( AL > ZERO) 

X PTR = X* T 
XPTI = Y*t 

CXPT = DC Mp LX (XPTR f XPTI ) 

CNUM = UCS1*CS1*C$1) + (CS2*CS2*C$1 ) ) *C DEX P( CX PT ) 

CS3 = ( (CS 1*CS1) - <C$2*C$ 2) ) 

CDENO* = CS3*CS3*C$3*CS3 
CRESLT = Z N UM /COE MO M 
U = DREAL(CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4,01 
yp - 2,1 

mmmmm = 22 

LCSC = 3 

SUM = 2.00568D+2 
SUMA = 2.01796D+2 
TIME = 0.55822D-1 
ANALYTICAL = 2.00568D+2 
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Case 48 




r n,rl \ n • 





SU8R0UT I NE VALUE (X ,Y .t , 1L,'J ,V , mmmmm j 
C THIS IS TARLE ENTRY NUMBER A8 OF LAPLACE TRANSFORMS 8Y SPIEGEL 

IMPLICIT REAL *3 ( At B * D-h,U-Z ) 

I MPLICIT C3 M P LEX* 1.6 (C) 



mmmmm 


s MMMMM f l 


ZERO = 


0 .D+0 


CS1 = 


DCMPtXCX, YJ 


XPTR = 


X* T 


xn = 


y*t 


CXPT = 


OCMPLX(X^TR,XPTI ) 


CNUM = 


COEXP ( CXP T ) 


CS4 = 


CS1*CSI *C$1 


CS5 = 


OCMP LX ( A L » ZE RO ) 


CS6 = 


CS5*C$5*C$5 


COENO M 


^ C$4 ♦ CSo 


CRESLT 


= :num/cdenom 



V = OIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X p = -4.44 

yp = 2.2 
MMMMM =23 
LOSC = 3 

SUM = 3.14683D+1 
SUMA = 3.14767D+1 
TIME = 0.48958D-1 
ANALYTICAL = 3.14683D+2 
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Case 49 




c at/ * f VJ at 
■^rjcos— + 



V3 



\^2 at 

8m — 2 — 



0 



- J«r/2 



} 



SU8R0U T I Nc VALUE (X , Y f T , AL ,U ♦ V , ) 

THIS IS TABLE ENTRY NUMBER 49 OF LAPLACE TRANSFORMS 
ILLICIT PEM_*8 ( A, 3, C-ht 3-Z ) 

I MPLICIT CO M PLEX* 16 i C) 

MMMMM = MMMMM v 1 
ZERO = O.O+O 
C SI = DC MPL X ( X f Y) 

XPTP = X*T 
XP T I = Y*T 

CXPT » CCMPLX(XPTR, XPTI) 

C S3 = CDEX^trx^T ) 

CNUM = CSl*CS3 
CS4 = C$1*CS1*CS1 
CSS = DCMPLXt ALtZEROJ 
CS6 = CS5*CS5*C55 
CDENOM = C S4 + CS6 
CRESL" = C NJ V/C OE N0 M 
U = 0 RE AL (C RE 5L T I 
V = DIMAG(CRESLT) 

RETURN 

END 



BY SPIEGEL 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
x F = .4.54 
y F - 2.2 
jMmmmm = 23 
Lose = 4 

SUM = 7.66825D+0 
SUMA = 7 . 92326 D+O 
TIME = 0.52416D-1 
ANALYTICAL = 7.66825D+0 
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Case 50 




K*' 



•* + 2e« f/a < 



>¥) 



SUBROUTINE V4LU£ ( X , Y , T ,4 L, U , V , 

THIS IS TA3LE ENTRY NUMBER 53 OF LA °L ACE 
IMdlICI T RC-AL*8 < A,8 ,D-H,0-Z> 

IMPLICIT COMPLEX* 16 (C) 

MM.MMM = mmmmm + 1 
ZERO = 0.0 + 0 
CSI = DCM°LX(X»Y) 

CS2 = CS1*CS1 
XPTR = X*T 
X FT I = Y* T 

CXPT = OCMPLXIXPTR.XPTI ) 

C S3 = CDEXP(CXPT) 

CWJM = CS2*CS3 
C S4 = C S1*C S1*C 31 
CS5 = OCMPLXI AL,ZERO) 

CS6 = CS5*CR5*CS5 
CDE.NOM = - S4 + C So 
CRESL T = CNUM/CDENOM 
U = 0 RE AL ( CRESL T ) 

V = DIMASICRESLT) 

RET J RN 

ENO 



TRANSFORMS 3Y SPIEGEL 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4.89 

y f = 2.3 

N1MMMM = 24 
LOSC = 5 

SUM = 8.34720D-1 
SUMA = 3.83818D+0 
TIME = O. 55796 D-I 
ANALYTICAL = 8 . 3472 OD-I 
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Case 51 





C-lat/2 







SUBROUTINE VALUE ( X, Y , T, 4 L , U, V , MM M mm) 

C THIS IS "ABLE ENT RY N'JM 0 £P 51 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLi: IT REAL*8 (AtBtO-H ,0-Z) 

ImolICIT CnMo LE X^16 ( C ) 

MMMMM =■ + 1 

ZERO = O.D+O 
CS1 = DCMPLX(X,Y) 

XPTR = X* T 
XPTI = Y*T 

CX D T = OC^RLXIXPTR, XPTI) 

CNUM = CDEXP(CXPT) 

CS4 = CS1*CS1*CSI 
CSS = D CMP LX (AL, ZERO) 

CS6 = CS5*CS5*CS5 
CDEMOM = CS4 - CSo 
CPESL* = C M'J M / COE NO M 
U = DREAL(CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

ENO 



AL = 1.25D-1 

T = 8 «D+0 

AA = 0,4 

= .4.44 
r 

y F • 2.2 

MMMMM =23 
LOSC = 3 

SUM = 3.25349D+1 
SUMA = 3.25*040+1 
TIME = 0.61204D-1 
ANALYTICAL = 3.25349D+I 
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Case 52. 




f 



~Za 




. y/3 at 




SUBROUTINE VALUE (X ,Y , T t AL,U , V , mmmmm ) 

C T"I S IS TABLE ENTRY NUMBER 52 OF LAPLACE TRANSFORMS BY SPIEGEL 

implicit real *s < a,3,o-h, o-zj 
I RPLICIT COMPLcX-16 (C) 

MMMMM = MMMMM + 1 
ZERO = 0.3+0 
CSX = OCMPLX(X,Y) 

XPTR = X*T 
X PT I = Y*T 

C X°T = DCMPLX (X^R, XPTI ) 

CS3 = COEXP(CX D T) 

CMJM = CS1*CS3 
C$4 = C$1*CS1*CS1 
CS5 = DCMP LX ( AL t Z ERG ) 

CS6 = CS5*C$5*CS5 
CDENOM -- CS4 - C 36 
CRESLT = C NIJ M /C OE MOM 
U = ORE AL (CRESLT ) 

V = OIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 



T = 8.D+0 



AA = 0.4 




MMMMM = 24 



LOSC = 4 



SUM = 8 . 33492D+0 
SUMA = 9.S6283D+0 
TIME = 0.67418D-I 



ANALYTICAL = 8 . 33^920+0 
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Case 53 



z 







eat + 2«” <l,/ * 



>fZat \ 



$ U953UTIME V*Ll)E ( X , Y , T t A L , U , V, MM mmmj 

T HI S IS TABLE EM'RY NUMBER 53 OF LA°LACE TRANSFORMS BY SPIEGEL 
I MPLI Cl T R FAL*8 ( A ,3 , D-H , G-Z ) 

I M PL ICIT COMPLEX* lt> (C ) 

MMMMM = MMMM.M + i 

ZEPO = 0,0 + 0 
CS1 = Dcr PLX(X,Y) 

CS2 = C$1*CS1 
XPTR = X*T 
X PT I = Y*T 

CXPT = DCMPLX(XP T R,XP^I) 




CDENOM = C$4 - CS6 
CRESLT = CNUM/CDSOM 
U = DR E AL ( CP ESLT ) 

V = DIMAG(CRESLT) 

R ETJRN 
END 



AL = 1.25D-1 



T = 8.D+0 



AA = 0,3 



X ? = -4.99 

y F = 2.3 



MMMMM = 24 



LCSC = 5 



SUM = 1.16806D+0 



SUMA = 2.25373D+0 



TIME = 0.70200D-1 



ANALYTICAL = 1.16806D+0 
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Case 5^- 




——(sin at cosh at — cosafsinhaf) 
4 a 3 



SUBROUTINE VALUE { X , Y r T , AL t U , V t mm mmm ) 

C THIS IS ""ABLE ETRY M'J v £ ER 54 OF L^IACE TRANSFORMS BY SPIEGEL 

IMPLICIT RF4L*8 ( A* B , 0-H , 3-Z) 

I V PL I C IT rnvo lex*1o (CJ 
MMMMM = MMM^M + L 
ZERO = 0.0+0 
CSX = DC*°LX(X,Y) 

X PTR » X* T 
XPTI = Y*T 

CXPT = CCMDLX (X°TR, xori) 

cnum = cdexpicxpt; 

CS4 = CS1*CS1*CS1*C$1 
CS5 = 0 CMP LX ( AL » ZERO! 

CS6 = 4.0+0*(CS5*C$5*CS5*CS5) 

COENOM = CS 4 + CS 6 
CRESLT = CNUM/COENOM 
U = OREAL(CRESLT) 

V * 0 I M AG ( CR ESL T J 

RETURN 

END 



A L = 1.25D-1 
T = 8 .D+0 
AA = 0,4 
X F = -4.44 

y F = 2 - 2 

MMMMM =23 
LOSC = 2 

SUM = 8.^92?2D+1 
SUMA * 8.U9272D+1 
TIME = 0.67732D-1 
ANALYTICAL = 8.49272D+1 
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Case 55 




-I 



r 4 + 4a* 



sin at sinh at 
2a 2 



SUBROUTINE VALUE ( X, Y, T, AL , U, V » MMMMm) 

C THIS IS TABLE ENTRY NUMBER 55 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REALMS ( A, 3 , 0-H ,3-ZI 
IM=>l ICIT r OVD L £ X ^lo ( C ) 

MMMMM = + L 

ZERO = 0,0 + 0 
CS1 = OC MO L X ( Xf Y ) 

XPTR = X**" 

XPTI = Y*t 

CXPT = DCM PLX(XP7R*XPTI ) 

C S3 = CDEX - ( CXPT J 
CNUM = C S1*C S3 
CS4 = CS1*C$1*CS1*CS I 
C$5 = DCMP LX ( AL > ZERO ) 

CS6 = 4.D+0*(CS5*CS5*CS5*CS5) 

CDENOM = CRA + C3o 
CRESLT « CN UM/CDt NQM 
U = OREAL(CRESLT) 

V = DI M AG (CRESLT i 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X F = -4.54 

y F = 2.2 

MMMMM = 23 
LOSC = 3 

SUM = 3.16447D+1 
SUMA = 3 . 16532D+1 
TIME = 0.59356D-1 
ANALYTICAL = 3.16447D+1 
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Case 56 




_L 

2a 



(sin at cosh at 



+ 



cos at sinh at) 



SUBROUTINE VALUE (X,Y,T, AL f U f V,MMMMM) 

C THIS IS TABLE ENTRY NUMBER 56 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT PEAL*8 ( A,B,D-H,0-Z) 

I HPL I C I T C0MPLEX*16 ( C) 

MMMMM = MMMMM > 1 
ZERO = 0.0+0 
CSI = DC MP LX ( X i Y ) 

CS2 = CSl-CSi 
X PTR = 

XPTI = Y*T 

CXPT = DCMPLX( XPTR, XPTI ) 

CS3 = COEXP(CXPT) 

CNUM = CS2*CS3 

CS4 = CSI*CS1*CS 1*0$ 1 

CS5 = DCMPLX ( AL i ZERO) 

C $6 = 4.D+C*(CS5*C$5*CS5*CS5) 

CDENOM = rs4 + CS6 
CPESLT = CNUM/ CDENOM 
U = OREAL(CRESLT) 

V = DIM AG( CR ESLT ) 

RETURN 

ENO 



AL = 1.25D-1 
T = S.D+O 
AA = 0.4 
X F = -4.89 
y F » 2.3 

MMMMM =24 
LOSC = 4 

SUM = 7.73369D+0 
SUMA = 9.23364D+0 
TIME = 0. 762C6D-1 
ANALYTICAL = 7.73369D+0 
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Case 57. 




cos at cosh at 



SUBROUTINE VALJE ( X t Y * T , AL r U» V , 

C this IS ~A6LE ENTRY NUMBER i>7 OF LA°LAC E TRANSFORMS BY SPIEGEL 

I HP LICIT REALMS ( A ,3 *0-H, O-Z) 

IMPLICIT C3MPLEX*16 (C) 



MMMMM = MMMMM + I 



ZERO = 


0.0+0 


CS1 = 


OCM°LX ( X * Y ) 


CS2 = 


CSI*CS1*CS1 


XPTR = 


X*T 


X PT I = 


Y 


CXPT = 


DC M P LX ( X PT R f X PT I ) 


C S3 = 


COEXP(CXPT) 


CNUM = 


CS2*CS3 


CS4 = 


C$1*CPI*CS1*CS1 


CS5 = 


DCMPLXI AL* ZERO) 


CS6 * 


^.Q+0*(CS5*CS5*C35 


CDENOM 


= CS4 + CS6 


CRESL * r 


= CNUM/CDENOM 


IJ = OP 


EAL(CRE3L X ) 


V = DI 


MAG(CRESLT) 


RETJRN 




END 





AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
Xp - -4,89 

y F - 2.3 

MMMMM = 24 
LOSC = 5 

SUM = 8.33730D-1 
SUMA = 3.35232D+0 
TIMS = 0.82576D-1 
ANALYTICAL = 8.3373CD-1 
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Case 58. 




^j(sinha« - sin at) 



c 



SUBROUTINE VALUE (X , Y , T .4 L, U , V, MMMmm) 

THIS IS T ABLE S' IT RY M U M 6 E 5 55 OF LAPLACE TRANSFORMS 8Y SPIEGEL 
IMPLICIT REALMS (Af5.C-H,0-ZJ 
IMPLICIT COMPLEX* 16 (C ) 



M V MMM s + i 



ZERO = 


G.0 + 0 




CS1 = 


OCMDLXIX, YJ 




X PTR = 


X*T 




XPTI = 


Y*T 




CXPT = 


DC m P L X ( XPTR 


f XPTI ) 


CNUM = 


CDEXP< CX or ) 




CS4 = 


CS1*C si*csi* 


CS1 


C$5 = 


OCMPLX* AL, ZE 


? n ) 


CS6 = 


CS5*CS5*CS5*CS5 


CDE'OM 


= C$4 - CS6 




C RES L^ 


= CNU M / CDENOM 



V = DIMAG(CRESLT) 

RETURN 

END 



AL = 1 . 25 D -1 
T = 8 .D +0 
AA = 0.0 
x F = - 0.04 
y p = 2.2 
MMMMM =23 
LOSC = 2 

SUM = 8 . 5 O 3 O 9 D+I 
SUMA = 8 . 5 - 30902+1 
TIME = O. 7 OO 90 D-I 
ANALYTICAL = 8 . 503 O 9 D+I 



150 



Case 59. 




S WOUT INE VALUE ( X, Y, T , AL, U, V, MMy mm) 

C THIS IS TABLE EMT^Y rjMBE^ 59 nr LAPLACE transforms BY SPIEGEL 
IMPLICIT REALMS ( A , 3 , J-h ,0-Z I 
IMPLICIT CD MO l EX* lb (C) 
mmmmm =. mmmmm + i 
ZERO = O.O+O 
CS1 = OCM^LXIX, Y) 

X pr R = X*T 
XPTI = Y*T 

C X D T = OCMPLX(XPTR , XPTI ) 

C S3 = CDEX p (CXPT j 
CNUM = CS 1*C $ 3 
CS4 = CS1*C$1*CS1*CS1 
CS5 = DCMPLXUL, ZERO) 

CS6 = C S5*CS5*C$5*C S5 
COENOM = CS 4 - CSfc 
CRESLT = CNUM /COE NOM 
J = DR EALICRESLTJ 
V = dimagicreslti 
RETURN 
END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
^ = -*. 5 * 
y F = 2.2 
MMMMM = 23 
LOSC = 3 

SUM = 3.20889D+1 
SUMA = 3 . 20897U+1 
TIME = 0.7S390D-1 
ANALYTICAL = 3.20889D+1 
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Case 60 . 




2^-{sinhat +• sin a*) 



c 



SUBROUTINE VALUE (X ,Y ,T f AL *U t V ♦ MMMmm ) 

THIS IS TABLE ENTRY NUMBER 60 jF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL *8 ( A f ) 

I MPL I C I T C0MPLcX*i6 (C) 

= MMMMM 4- I 



ZERO = 0.0+0 
CSI = DCMPLXI X, YJ 
CS2 = CSi*CSl 
X PTR = X*T 
XPTI = Y*T 

CXPT = DCMPLX(XPTR,XPTI ) 

CS3 = C DEX P< CXPT J 

CNLM = C S2* C S3 

CS4 = CSl* r $l*CS 1*CS I 

CS5 = DCM? LX( AL,ZERD> 

CS6 = C S5*C S5*C$5*C$5 
C OENHM = CS 4 - CS 6 
CRESLT « CNUM/CDENOM 
U = D°EAL(C D ESLT) 

V = DIMAGICRESL7) 

RETURN 

END 



AL = 1.25D-1 
T = 8,D+0 
AA = 0.3 
X = -4,54 
y = 2,2 
MMMMM = 23 
LOSC = 4 

SUM = 8 . 06669D+0 
SUMA = 8.32767D+0 
TIME = 0 . 87672 D -1 
ANALYTICAL = 8 . 06669 D+O 
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Case 61 



^(cosh at + cos at) 



SUBOUT INE VALUE ( X , Y , T , A L , U , V f MMMMm) 

C THIS IS TABLE EN r RY NUMBER 61 OF LA °L AC E TRANSFORMS BY SPIEGEL 

IMOLi:iT REAL*8 (A,8fO-H f O-Z) 

IMPLICIT r "lMPLcX*l6 (C ) 

MMMMM = MMMMM + 1 
ZERO = 0,0+0 
CS1 = DC MO L X ( X t Y ) 

CS2 = CSl*C$i*C$i 
XPTR = X*T 

XPTI = Y*t 

CXPT = DCMPLX(XPTRfXP T I ) 

C$3 = CDEXP(CXPT) 

CNJM = CS2*CS3 

C$4 = CSL*CSi*C$I*CSL 

CS5 = 0 C M d l X ( AL » Z ER3 ) 

CSS = CS5*CS5*CS5 *C$5 
CDEN3M = CS4 - C S 6 
CPESL7 = CNUM/COENOM 
U = DPEAL ( CP ESL X ) 

V = DIMAG(CRESLT) 

RETURN 

END 




AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
X F = -4.99 
y F = 2.3 
MMMMM =24 
LOSC = 5 

SUM = 1.04169D+0 
SUMA = 2.13326D+0 
TIME = 0.93314D-1 
ANALYTICAL = 1.04169D+0 
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Case 62 




€ -bt - 

2(6 - <x)yfr& 



SUBROUTINE VALUE (X, Y, T, A L , BE, U, V, MMMMM, NCALL) 

C THIS IS "ABLE E k r*Y NUMBS® 62 OP LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT R £AL*8 ( A , 8 , D-H ,0-Z) 

IMolICIT COMPLEX* 16 (C) 

DIMENSION NO ALL 15) 

MMMMM = MMMMM + 1 

ZERO = O^O+O 
CSI = DCMPLX(X,Y) 

CS2 = DC MP LX< AL , ZERO) 

CS3 = DCMPLX (3E,ZER0) 

X PTR = X* T 
XPTI = Y*r 

CX °T = DCMOLX {XPTRtXPTD 
CNUM * CDEXP(CXPT) 

CS4 = CSI + C$2 
CS5 = CSI + CS3 

CALL CPOWER {CS4,CS6,5.0-1 ,1 , NCALL) 

CALL roDWE®. <CS5,CS7, 5,0-1, 2, NCALL) 

CDENOM = CS6 + CS7 
CRESLT = C MUM /COE NOM 
U = DPEAL(CPESL") 

V = DIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
BE = -1.25D-1 
T = 8 .EH-0 
AA = 0,3 
X F = .4,52 

y F = 2.2 

MMMMM = 23 

Lose = 5 

SUM = 1.17209D-1 
SUMA = 6 . 50880D-1 
TIME = 0.13301D+0 
ANALYTICAL = 1.17209D-1 
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Case 63* 



1 

» + a 



erf y/at 

V* 




SUBROUTINE VALUE (X, Y, T, Al, U, V, MMMM* f NCA LL) 
THIS IS T ABLE ENTRY NUMBER 63 OF LAPLACE TRANSFORMS 
IMPLICIT C EAL*8 ( A,B f D-H.Q-Z) 

ILLICIT COMPLEX* 16 (C) 

DIMENSION NCALLC5J 

+ I 

ZERO = O-D+O 
C SI * DCMOLXI X,Y) 

C$2 = DC MP LX ( AL » ZERO ) 

XRTR = X* T 
XPTI = Y*T 

CXPT = DCMPLX(X p ^R, XPTI ) 

C NUM = CDE XP(CXPT) 

C S3 = CS1 + CS2 

CALL CPQWER ( CS 3, CS4, 5 .C- 1, 1, NCA LL) 

CDENOM = CS1*CS4 
CRESLT = CN U M /CDENOM 
U = ORE AL ( CRESLT ) 

V = D I MAG (CRESLT) 

RETURN 

END 



BY SPIEGEL 



AL = 1.25D-1 
T = 8 .D+0 
AA = 0.2 
X p = -4.64 

y F * 2.2 

MMMMM =23 
LCSC = 5 

SUM = 2.38352D+0 
SUMA = 2.6473^0+0 
TIME = 0.85202D-1 
ANALYTICAL = 2.38352D+0 
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Case 64 



e at er/ '/at 

\Tc i 




SUB ROUT I ME VALUE (X, Y, T, AL, U 
THIS IS r A 6L t ENTRY NUMBER 64 OP 
IMPLICIT REALMS ( A, 3, C-h, 0-Z) 
IMPLICIT COMPLEX^ 16 (CJ 
0 IMENS ICN NC ALL (5 ) 

MMMMM = MMMMm + 1 
ZERO = 0.0+0 
CS1 = DC MPLX ( X , Y ) 

C$2 = DCMPLXCAL, ZERO) 

X D ‘R = X*T 
X pr I — Y*T 

CXPT = DCMOLX(XPTR,X£>TI ) 

CNUM = CDEX?( CXPT ) 

CALL CPCWER (CSl»CS3 f 5.0-lfIfNCALL) 
CS4 = CS1 - CS2 
COEMOM = CS3*C$4 
CRE SL T = C MUM/ COE MOM 
U = ORE AL (CRE SL T) 

V = DIM AG( C° ESLT ) 

RETURN 

END 



V, MM«MM r NCALL) 

LAPLACE TRANSFORMS BY SPIEGEL 



AL = 1.253-1 

T = 8.D+0 

AA = 0.4 

X = -4.44 
F 

y F - 2.2 

MMMMM =23 
LOSC = 5 

SUM = 6.47907D+0 
SUMA = 8.56802D+0 
TIME = 0.97552D-1 
ANALYTICAL = 6.47907D+0 
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Case 65. 









SUBROUTINE VALUE (X, Y, T, AL t 8E» U, V, MMMMM, NCALL) 

THIS IS T A 3LE ENTRY NUMBER 65 OF LAPLACE TRANSFORMS 8Y SPIEGEL 
IMPLICIT REAL*8 < A, 3 ? C-h,3-Z) 

IMPLICIT c nMD LEX* 16 f C ) 

DIMENSION NC ALL (5 ) 

MMMMM = MMMMM + L 
ZERO = O.O+O 
CS1 = DC M»LX ( X » Y J 
CS2 = DC MP LX( AL f ZE RO) 

CS3 * DCMPLXI BE, ZERO! 

XF T R = X*t 
XPTI = Y*T 

CXPT = CCMPLX(XPTR ,XPTI ) 

CNUM = CDEXP(CXPT) 

C S4 = CS1 - CS2 

CALL rpnwER ( CS 4, CS 5, 5- 0- 1, 1, NC ALL ) 

CDENOM = CS 3 + CS 5 
CRESLT = CNUM/COENOM 
U = DREAL ( CRESLT J 
V = DI MAG (CRESL7J 
RET N 
END 



AL = 1.25D-1 

BE = 3.535533906D-1 

T = 8.D+0 

AA = 0.3 

X p = -*.5* 

y F = 2.2 

MMMMM = 23 
LOSC = 5 

SUM = 1.31286D-1 
SUMA = 9.058^90-1 
TIME = 9.0116D-2 
ANALYTICAL = 1.31286D-1 
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Case 66 



i 




Jq («0 



SUBROUT INE VALUE (X, Y, T, AL, U, V, M MMHM, NCALL) 

C THIS IS TABLE E IT RY NIJM 3 £P. 66 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL*8 ( A , 3 , 0-H ,0-Z) 

IMPLICIT COMP L EX* 16 l C ) 

01 MENS I ON NCALL (5 i 
MMMMM = MMMMM > 1 

ZERO = 0,0+0 
CS1 = DCMPLX(XtY) 

CS2 = DCMPLX(ALtZERO) 

XP”R = X* T 
X PTI = Y*“ 

CXPT = OCMPLXI XPTR, XPTI ) 

CNUM = COE X D ( CX PT ) 

C S3 = (CS1*CSI) + (CS2-CS2) 

CALL CPOWEP (CS3 f C0ENGM ? 5.D-i, i i NCALL) 

CRESLT = CNUM/COENOM 
U = DP E AL( CRESLT) 

V = DIM AG( CRESLT ) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.2 
Xp = -k.6U 

y F = 2.2 

MMMMM =23 
LOSC = 5 

SUM = 7.65198D-1 
SUMA = 1.22916D+0 
TIME = 0. 91256D-1 
ANALYTICAL = 7 . 65198 D-I 
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Case 67 



1 



M* 0 




SUSROU 
THIS I 
I MPLIC 
I VPLIC 
0 1 M E'I S 
MMMMM 
Z £R 3 = 
CS1 = 
CS2 = 
XPTR = 
XPTI = 
CXPT = 
CNUM = 
CS3 = 
CALL C 
CRESLT 
U = 3® 
V = DI 
RE T URN 
ENO 



INE VALUE (X, Y, T, AL, Ui 
S TABLE ENTRY NUMBER 67 OF 
IT P EAL-a ( A, 3, C-H, O-Z) 

I T C2^°LEX- 16 ( C) 

ION MC ALL ( 5 ) 

- yyyMy + 1 

0*D + 0 

OCMPLXI X, Y ) 

OCMPLXI AL, ZERO) 

X*T 

Y*T 

OCMOLX (XPT^ ,XPT I J 
COE XP(CXPT) 

( C$l*CSl) - (CS2*CS2) 

POWER (CS3, CCENOM, 5.0-1, 1,NCALL) 
= CNUM/C3EN0M 
E AL (CRESLT) 

MAG(CRESLT) 



V, MMMMM, NCALL) 

LAPLACE T R AN SFORMS BY SPIEGEL 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.3 
x F = -*.99 

y F = 2.3 

MMMMM = 2* 

LOSC = 5 

SUM = 1.26607D+0 
SUMA = 2.31C76D+0 
TIME = 0 *102910+0 
ANALYTICAL = 1.26607D+0 
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ase 68 



n > - 1 



a*J H (at) 




C 



SUBROUTINE VALUE (X, V, T, AL » NPOWER, U, V* MMmmm. 
THIS IS TABLE ENTRY NUMRER 68 np LAPLACE T RANS F0R w S 
IMPLICI T R E A L * 8 ( A, 3 , 0-H f >Z ) 

I ELICIT C0MPLEX«16 ( C) 

DIMENSION NC ALL ( 5 ) 

MMMMM = MMMMM + 1 
ZERO = 0,0+0 
C SI = DC MP L X( X * Y ) 

CS2 = DCMP LX ( AL, ZERO) 

XPTR = X* T 
XPTI = Y*T 

CXPT = CCMOLX( XPTR, XPTI ) 

C S3 = (CS1*CS1) + (C$2*CS2) 

CALL CPOViER (CS3»CDENC M t5.D-Ifl #NCALL) 

C$4 = ( C DEMON - CS1)**NP'"HEP 

CNU^ = CSA*C0EXO( CX^T) 

C RE SL T = CNU M /COt NOM 
U = ORE AL ( C RESLT ) 

V = DI MAGICRESLTJ 

RETURN 

END 



NCALL) 

BY SPIEGEL 



AL - 1.250-1 
N POWER = 3 
T = 8.D+0 
AA = 0.4 
X F = -2.16 

y F = i.6 

MMMMM =17 
LOSC = 1 

SUM = 3 . S2096D-5 
SUMA = 3.82096D-5 
TIME = 0.8122UD-1 
ANALYTICAL = 3.32090D-1 
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Case 6S. 



£ 




o" I n (at) 



SUBROUTINE VALUE (X, Y» T, Al , N^PWER » U, V, NCALL) 

THIS IS TABLE ENTRY NUNS ER 69 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL-8 < A , 3 , O-H ,Q-Z > 

I ^PLICIT CC'^PL EX* lo ( CJ 

DIMENSI ON NCALL (5) 

= M'-'MMM + 1 

ZERO = 0 .0+0 
C SI = OCMPLXI X,Y) 

CS2 = DCMOLXI AL,ZERJ) 

X PTR = X*^ 

XPTI = Y-T 

CXPT = CCMPLX (X 3 TRt X DT I ) 

C S3 = (CSl*C$l) - ( C5 2*C5 2) 

CALL CPOUER (CSStCDENfM, 5.0-1 1 1 jNCALL) 

CS4 = (CS1 - CDLNOM)**npowE* 

CNUM = CC4*C0EXP{ CX°T ) 

CRESLT = CMUM/C06NJM 
U = DR E AL ( C RES LT J 
V * 01 MAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 



NFCWER - 3 



T = S.D+O 



AA = 0.4 




MMMMM = 1? 



LCSC = 1 



SUM = 4.32979D-5 
SUM A = 4 . 329790-5 
TIME = 0.74672D-1 



ANALYTICAL = 4.329770-5 
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Case 70. 



,t>is — 



Jo(aVT(t + 2b)) 




“1 



SUBROUTINE VALUE (X, Y, T, At » 3E, U, V, 
THIS IS TABLE ENT^Y NUMBER TO Of LAPLACE 
IMPLICIT * E*L*8 ( A ,3 ,0-H , O-Z ) 

IMPLICIT C0MPLEX*16 (C) 

DIMENSION NC ALL t 5 J 
M MM MM = MMMMM + I 

ZERO = 0,0+0 
CS1 = DCMPLX(X,Y) 

CS2 = D CMP LX ( AL ,ZERO ) 

CS3 = DC M PLX < 8£ » ZERO ) 

XPTR = X*T 
X PT I = Y*T 

CXPT = DCMPLX IX P7 R»X OT I ) 

CS4 = (CS1-CS1) + ( C S 2*C S2) 

CALL CPOWEP (CS4, COENOM, 5.D-1, I,NCALL) 
CS5 = (CSl - COE-OM)*CS3 
CNUM = C0EXP(CS5)-CDEX?(CX^ T ) 

CRESLT = CNMM/CDfcNOM 
U = DREAL(CRESLT) 

V = DIM AG( CR ESL T ) 

RETURN 

END 



MMMMM. NCALL) 
TRANSFORMS BY 



SPIEGEL 



AX = 1.25D-1 
BE = -1.25D-1 
T = 8 . D+0 
AA = 0.2 




MMMMM =23 

LCSC = 5 



* SUM = 7.72088D-1 
SUMA = 1.23U85D+0 
TIME = 0.8291^D-1 
ANALYTICAL = 7.72088D-1 
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Case 71 



t > b 
t < b 



„?•■[*] = 1 



— J J 0 (a\/l^~—b^) 
0 



SUBROUT INE VALUE (X, Y, T t <L, RE, 'J, V, MM MM M f NCALL ) 

THIS IS TABLE E'JTRY NUMBER 71 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL *8 ( A, 3,0-H»3-Z) 

IMPLICIT C0 M PLEX*io (C) 

0 IMENSI ON NC ALU 5) 

MMMMM = MMVMM + 1 
ZERO = 0,0+0 
CSI = 0 C M P L X ( X, Y) 

CS2 = D CMP LX ( AL t ZERO) 

C S3 = DC MPL X( -BE » ZERO) 

XPTP = X*T 
X PTI = Y* T 

CXPT = DC M° LX l XPTRtXPTI ) 

CS4 = ( CS i *C S i ) + <C$2*CS2) 

CALL CPOWER <CS4,CDENCM # 5.D-It 1, NCALL) 

CNUM = CDEXP(CDEN0M*CS3)*CDEXP< CXP T ) 

CRESL" = CNUM/CDENOM 
U = DRE AL(CRESLT) 

V = DIMAG(C-ESLT) 

RE T J RN 
END 



AL = 1.25D-1 
BE = -1.25D-1 
T = 8.D+0 
AA = 0.2 
X F = -U.6k 

y F = 2.2 

MMMMM =23 

Lose = 5 

SUM = 7.65252D-1 
SUMA = 1.2U690D+0 
TIME = 0.86658D-1 
ANALYTICAL = 7.65252D-1 
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se 72. 



1 

(*- -h a 2 ) 3/2 







a 



SUBROUTINE VALUE (X, Y, T, , U, V, MMMMM, NCALL) 

C THIS IS TABLE entry NU^RER 72 OF LAPLACE TRANSFORMS 8Y SPIEGEL 

IMPLICIT RriL^a ( A, 9 ,0-H, 3-Z ) 

I M PL I C I T COMPLEX* 16 < C) 

DIMENSION NC A LL ( 5) 

MMMNM = MMMMM > 1 

ZERO * O.D+O 
XPTR = X *T 
X P' r I = Y*T 

CXPT = DCMPLX(XP T R,XPTI) 

CNUM = ODE >P( CXP T ) 

X REAL = <X*X) - <Y*Y) + <AL*AL) 

XI MAG = 2.D + 0*X*Y 

CS1 = DCM°L X ( XR cA L» X I MA G ) 

CS2 = CSI*C$I*CS1 

CALL CPQWER ( CS2 , COENCV, 5.0-1, 1, NCALL ) 

CRESLT = CUU M /CQENOM 
U = DREAL(CRESLT) 

V = 0 1 MAGI CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
Xp = -k.M 

y F = 2.2 

MMMMM =23 
LOSC = 3 

SUM = 2 . 81632D+1 
SUMA = 2. a 1715D+l 
TIME = O.IIU 76 D+O 
ANALYTICAL = 2.81632D+1 
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( 8 2 + fl 2)3/2" 



tJ Q (at) 




SUBROUTINE VALUE (X, Y, T f A L , (1, V, MMM M M f MCML) 

C THIS IS TABLE ENTRY N U M BER 73 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REALMS (A,BfC-h»0-Z) 

IMPLICIT COMPLEX- 16 ( C) 

□ HEMS ION N C ALL ( 5) 

M MMMM = MMMMM + i 
ZERO = C.D+O 
CS1 = DCMPLX(XtY) 

XPTR = X* T 
XPTI = Y*T 

CXPT = DCMPLX (X°T R t X D T I ) 

CS2 = COEXP(CXPT) 

CNUM = CS1-CS2 

X REAL = (X -X ) - (Y*Y) + ( AL *A L ) 

XI MAG = 2.0+0*X*Y 

CS3 = DC MO LX(XREALtXIMAG) 

CS4 = CS3*CS3*CS3 

CALL CPOV.ER (C54, COENCM, 5.0-1, 1, NCALL) 

CRESLT = CNUM/COENOM 
U = DREAL(CRESLT) 

V = OIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.4 
X p = -4.89 
y F = 2.3 

MMMMM =24 
LOSC = 4 

SUM = 6.12158D+0 
SUMA = 8.3^9980+0 
TIME = 0.12620D+0 
ANALYTICAL = 6.121580+0 
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Case 7k 



(<*0 ** atj^at) 




SUBROUTINE VALUE (X, Y, T f Al, U t V, MMMMM, NCALL) 

C THIS IS TABLE E JTRY NUMBER 74 OF LAPLACE TFANSFQRMS BY SPIEGEL 

IMPLICIT REAL*3 ( A, 3, D-h,D-Z ) 

IMPLICIT COMPLEX* 16 (C) 

DIM5NS ION NC*LL( 5) 

MMMMM = MMMMM + 1 
ZERO = O.O+O 
CS1 = DC M D LX { X » Y ) 

XPTR = X*~ 

X PT I = Y*T 

CXPT = CCMOLX (XPTRtXPTI) 

C 52 = CDEXP(CXPT) 

C S3 - CS 1*CS1 
CMJM = C52 *CS 3 

XREAL = ( X*X) - ( Y*M) + ( AL* AL ) 

XIMAG = 2. 0+ 0*X * Y 

CS4 = DCMPLXIXREALtX IMAG) 

C S5 = C S4*C S4*C S4 

CALL CPOWER ( C S5 ♦ COENOM, 5.0-i, 1, NCALL ) 

CRESL* = CNUM/CDENOM 
U = 0REAL(CRE5LTJ 
V = DIM AG( CRESLT) 

RETURN 

ENO 



AL = 1.25D-1 
T = 8 .D+0 
AA =0.4 
X p = -4.89 

y F - 2 -3 

MMMMM = 24 
LOSC = 5 

SUM = 3.251^7D-1 
SUMA = 3.72547D+0 
TIME = 0 . 10772D+0 
ANALYTICAL = 3.251 U 7D-1 
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Case 75 



1 

(,2 _ „ 2 ) 3/2 



(/,(0i) 




a 



SUBROUTINE VALUE (X, Y, T, AL, 'J, V, MMMMM, NCALL) 

C THIS IS T ABL E E JT RY NUMBER 75 '•F LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REALMS ( A, 0-2) 

IMPLICIT COMPLEX* 16 (C) 

01 MENS I ON NCALL (5 ) 

MMMMM = M M M MM + I 
ZERO = 0.0+0 
XPTR = X*" 

XPTI = Y* T 

CX°t = DC.mplX(X?TR, X°T I ) 

C NUM = COE XPICXPT ) 

XPEAL = ( X*X) - <Y*Y) - (AL*AL ) 

X I MAG = 2 . D+0«X *Y 

CS1 = DCMPL X( XREALtXIMAG) 

CS2 = C$1*C$1*CS1 

CALL CPOWER (CS2, COENOM, 5.0-1, 1, NCALL) 

CPE SL T = CMUM/COE.NOM 
•J = 05 EAL(CPESLT) 

V = OIMAG(CRESL’) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.5 
X F = -^.3^ 

y F = 2 ' 2 
MMMMM =23 

LOSC = 3 

SUM = 3.617C2D+1 
SUMA = 3 . 638 A 6 D+I 
TIME = 0.11C89D+0 
ANALYTICAL = 3.61702D+1 
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Case 76. 




tl^at) 



SUB^CUJT INE VAL'JE (X, Y, T, AL , U, V, MMMMM, NCALL) 

THIS IS ’’’ABLE ENTRY NUMBER 76 OF LA°LACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL*8 ( A, 3 , O-H, O-Z ) 

IMPLICIT C3M?LEX*16 tC) 

0 IM6NSI ON NCALL (5; 

+ 1 

ZERO = O.D+O 
CS1 = DCMOLXf X, Y) 

XPTR = X*T 
XPTI = Y*T 

CXPT = DC^LX< XPTR, XPTI ) 

CS2 = CDEX c (CX p T ) 

CNUM = CS 1* CS 2 
XRE*L = (X*X) - ( Y*Y) 

XI*AG = 2 * D +0 *X *Y 
C S3 = OCMPLX (XREAL ? X I M A 6 ) 

C$4 = CS3*CS3*C$3 
CALL C POWE 0 (CS4, COSNO M , 

CRESLT = CNU m/CDENOM 
U = DP E AL ( C c ESL~ ) 

V = OIMAG(CRESLT) 

RETURN 
END 



i 4L*AL) 



5.0-lt ly NCALL) 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.5 
X p = -*.79 
y F - 2.3 

MMMMM =24 
L05C = 4 

SUM = 1.012S5D+1 
SUMA = 1.52218D+1 
TIME = 0.13250D+0 
ANALYTICAL = 1.01285D+1 
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Case 77 • 



7 0 (oO + atltiat) 




SUBROUTINE VALUE (X, Y, T, AL t U, V, MMMMM, NC ALL ) 

C THIS IS TABLE ENTRY NUMBER 77 OF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT 3 E AL *8 ( A, d, 0-F,0-Z J 

I MPLI C I T COUPLE X* 1 6 (C) 

. DIMENSION NC ALL ( 5 J 
MMMMM = mmmMM ♦ 1 
ZERO = 0,0+0 
CSl = dcmdlxi X t YJ 
x P-R = X* T 
XPTI = Y*T 

CXPT = CC M P LX(XPTP fXPTI ) 

CS2 = CDEXP(CXPT) 

C S3 = CS1*CSI 
CNUM = CS2 * r S 3 

XREAL = <X*X) - ( Y *Y ) - <AL*AL) 

XIMAG = 2.Q+0*X*Y 

CS4 = DCMPLX(XREAL, XIMAG) 

CS5 = CS4*CS4*CS4 

CALL CPOWEP (CS5, C DENOM , 5.D-1, I, NCALL) 

C RES L T = CNUM/ CDEN3 M 
U = DRE A L( C °E SL T) 

V = DIM AG( C° .ESL T i 
RETURN 
E tO 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.5 
Xp = -4.79 

y F = 2-3 

MMMMM = 24 
LOSC = 5 

SUM = 1.93122D+0 
SUMA = 7.71705D+0 
TIME = 0.11S07D+0 
ANALYTICAL = 1.83122D+0 
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Case 78 



coa 2 '/at 




SUBROUTINE VALUE (X f Y, T f AL , U, V, MMMMM. NCALU 
C THIS IS TABLE ENTRY NUMBER 81 OF LA°LACE TRANSFORMS BY SPIEGEL 

IMPLICIT R EAL*8 ( A , R , 0- H , D-Z ) 

IMPLICIT COMPLEX *16 <C) 

DIMENSION NC ALL ( 5 J 
MMMMM = MMMMM + X 
ZERO * 0.D+0 
C SI = 0 C M d L X ( X, Y) 

CS2 = DCMP LX ( AL ,ZERO ) 

CS3 = -C$2 /CS1 
XPTR = X*T 
XPTI = Y*T 

CXPT = DCM^LX (XPTR.X ot IJ 
CNUM = CDEXP(CSj) -CDEXP(CXPT) 

CALL CPCWER CCSlt CDENOM, 5-D-lt 1»NCALLJ 
C R£SL T = CNUM/CDENOM 
U = DREAL(CRESLT) 

V = D I M AG ( CRES LT ) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.2 
Xy = -4.64 

y F = 2 - 2 

MMMMM =23 
LOSC = 5 

SUM = - 8 . 30093D-2 
SUMA = 4.47194D-1 
TIME = 0.96148D-1 
ANALYTICAL = -8.30093D-2 
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Case 79 




f — a/t 




c 




MMMM ‘4 = MMMMM 4* 1 
ZER3 = 0 . D+0 
CS1 = DCMPLX(X,Y) 

CS2 = OCMO LX(AL ,ZERJ) 

C$3 = -CS2/CS1 
XPTR = X*T 
XPTI = Y*T 

CXPT = 0C MO LX (X PT R ? X°T I ) 

CNUM = COEXF(CS3J ^COEXP{ CXP^ J 
CS4 = CS1*CS1*0S1 

CALL C POWER (CS4? C CENOM , 5 .0- 1 1 1 » NCALL ) 
CRESLT = CNUM/CDENOM 
U = DPEAL (C°E$LT) 

V =r OIMAG(CRESL T J 

RETURN 

END 



AL = 1.25D-1 



T * 8. D+O 



AA = 0.2 




MMMMM =23 



LOSC = 5 



SUM = 1.45103D+0 
SUMA = 1.82676D+0 



TIME = 0.10013D+0 



ANALYTICAL = 1.4<103D+0 
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Case 80 



■1 






n > -1 



£ 



(!)■'%. (2^) 



SUBPOUTINE VALUE (X, Y, T, AL, N»rWER, U, V, NCALL ) 

C THIS IS TABLE ENTRY NUMBER 83 OF LAFL4CE TRANSFORMS BY SPIEGEL 

I^PLIC IT o EAL*8 ( A,3,0-H,0-Z) 

IMPLICIT CO ydlEX * 16 (C) 

OIMENSION NCALL (5) 

MMMMM = MNM«M + i 

ZERO = 0.0+0 
CS1 = OCNPLXI X, Y) 

CS2 = OCMPLX (AL , ZERO ) 

XPTR = X*T 
X°TI = Y* T 

CXPT = OCMPLX (X’TR.XP* I ) 

CNUM = COEXP( CS3) *CDEXP(CXP-) 

COENO* = C S 1 **( NPO WE° + 1) 

CBESL T = CNUM/C3EN0N 
U = OPEAL(CRESLT) 

V = DIM AG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
N POWER = 3 
T = 8.D+0 
AA ' 0.3 
Xp = -4.11 
yp = 2.1 
MMMMM =22 
LOSC = 4 

SUM = 6.60189D+1 
SUMA = 6.60189D+1 
TINE = 6.9758D-1 
ANALYTICAL = 6.60139D+1 
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Case 81 . 




SUB ROUT I ME VALUE (X, Y, T, i,|_, U, V, MMMMM, NCALL) 

C THIS IS T4RLE ENTRY NUMBER 3A OF LA°LACE TRANSFORMS BY SPIEGEL 

IMPLICIT 0 cAL *3 < A,8 ,D-rl,D-Z ) 

I RPL I C I T CO MO LEX* 16 ( CJ 

DIMENSION NC A LL( 5 ) 

MMMMM = MMMMM I 
ZERO = O.O+O 
CSI = DC M°L X( X » Y) 

CS2 = OCMPLXIAL.ZERDJ 
XPTR = X* T 
XPTI = Y*T 

CXPT = OCMPLX (XPTR.XPT I ) 

CALL CPOWFR (CSI ,CCENC U , 5.0-1 ,1 ,NCALL) 

: S3 = -CS2*CDEN3M 

CNUM = CDEX P( CS3 ) *COEX-p (CXPT ) 

CRESLT = C MUM/COc NON 
U = DP EAL(CRESLT) 

V = D I M AG ( CR ESL T I 

RETURN 

END 



AL = 1.25D-1 
T * 8.D+0 
AA = 0.2 
X F = -U.6U 
y F = 2.2 

MMMMM =23 
LOSC = 5 

SUM = 1.9937^D-1 
SUMA = 6.27^dD-l 
TIME = 0.93^700-1 
ANALYTICAL = 1.99374D-1 
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Case 82. 




- 2 — ,-«■/<! 
2v'rl 3 



SUBROUTINE VALJE (X, Y, T, AL , U, V, MMMMM, NCALL) 

C THIS IS T ASLE ENTRY NUMBER 36 n F LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT REAL*8 < A,B ,D-H,3-Z) 

IMPLICIT CO«PLEX*l6 (CJ 
DIMENSION NCALL ( 5 ) 

MMMMM = MMMMM + 1 
ZERO = O.D+O 
C SI = DCMPLXI X,Y) 

CS2 = DC MO LX(4LtZ5R0) 

XPTR = X* T 
XPTI = Y*T 

CXPT = CCMPLXIXOTR, XPTI ) 

CALL C POWER (CS1, C33.5.C-1.1 , NCALL) 

CSA = - CS2*C S3 

C RES L T = CDcX P ( CS A )*CDEXP( CXPT ) 

U = DREAL(CRESLT) 

V = DIM AG(CRESLT) 

RET'JRN 
E NO 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.2 
X p = -5-09 

y F = 2.3 

MMMMM = 24 
LOSC = 6 

SUM = 1.55?6lD-3 
SUMA = 3.03773D-1 
TIME = 0.10197D+0 
ANALYTICAL = 1.55761D-3 
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erf (a/2 \/7 ) 




SUBROUTINE VALUE (X, Y, T, AL, U f V, MMMMM, NCALL) 

THIS IS ^ABLE ENTRY NUMBER 66 OF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REALMS ( A , B , 0-h, O-Z ) 

IMPLICI T COMPLEX- 16 ( C) 

0 I MENS ION NCALL (5) 

MMMMM = MMMMM + 1 

ZERO = O.O + O 
ONE = 1.0+0 
CS1 = DC wp LX( X, Y) 

CS2 = 0CMO{_X< AL ,ZERO) 

CS3 = DCMPLX (ONE, ZERO ) 

X PTR = X- T 
XPTI = Y*T 

CXPT = OCMPLX (X°TR,XP T I) 

CALL CP OWE R (C$1 ,CS4, 5.0-1 ,1 , NCALL) 

CNUM = ( CS 3 - COEXP ( -C S2*C S4) ) -COE XP (CXPT) 

CDENQM = CSl 
CRESLT = CNUM /COE NON 
U = DR EAL (CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

ENO 



AL = 1.25D-1 
T = 8 . D+0 
AA = 0.2 
Xp = -4.6* 

y p =2.2 

MMMMM = 23 
LCSC = 5 

SUM = 2.49298D-2 
SUMA = 8.07U80D-2 
TIME = 0.93158D-1 
ANALYTICAL = 2.49292D -2 



Case 84. 




erfc (0/2VT) 



SUBROUTINE VALUE IX, Y, T, AL, U, V, MMMMM, NCALL) 

THIS IS ’"ABLE ENTRY NUMBER 37 "> F LA D LACE TRANSFORMS BY SPIEGEL 
IMPLICIT R EAL* 3 I A ,3 , O-rl , 0-2 ) 

IMPLICIT COMPLEXES I C ) 

DIMENSION NCALLI 5 I 
MMMMM = M^MMM * 1 

ZERO = 0.0 + 0 
CSI = DCMPLX(X,Y) 

C $2 = DCMDLXI AL.ZcRO) 

XPTR = X*~ 

XPTI = y* t 

CXPT = OCMPLXIXPTR.XPTI) 

CALL CPOWER (CSl, CS4, 5.C-1, l.NCALL) 

CMJM = CDEXP{-CS2*CS4)*CD£XP(CXPT) 

CDENOM = rsi 
C RES L T = C M UM/ COENO M 
U = OREAL(CRESLT) 

V = DIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T = 8.D+0 
AA = 0.2 
X p * -4.64 

y F = 2.2 

MMMMM = 23 
LOSC = 5 

SUM = 9.75070D-1 
SUMA = 1.35850D+0 
TIME = 0.99658D-1 
ANALYTICAL = 9.75071D-1 
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Case 85 




,»<»(+«> tr { c / 4. -2—^| 

V 2 -ft) 



SUBROIJT INE VALUE (X, Y, t. AL > RE, a, V, 
THIS IS TABLE ENTRY NUM3ER 88 OF LAPLACE 
IMPLICIT REAL *8 < A, 3, C-h.O-Z ) 

I MPLICI T COMP LEX* 1 6 ( C J 
DIMENSION N CALL I 5 ) 

MMMMM = MMMMM + 1 

ZERO = 0. 0+0 
CS1 = DCMP L X( X, Y) 

CS2 = DCMPLX 1 AL, ZERO) 

XPTR = X* T 
XPTI = Y*T 

CXPT = DCMOLX(X°7R,X = Tn 
C S3 = DCMDLXI Be .ZERO) 

CALL C POWER ICS 1,03 A, 5.0- 1, L,MC ALL) 

CNUM = CDEX P I-CS2 *C.S 4 ) *CDEXP ( CXPT) 

CS5 = C S3 + C$4 
C0EN3M = CS4*CS 5 
C RE SL T = CNIJM/COENOM 
U = DREAL(CPFSLT) 

V = DIMAG(CPESLT) 

RETURN 

ENO 



MMMMM, NCALL ) 
TRANSFORMS BY SPIEGEL 



AL = 1.25D-1 

BE = 3.4574089063-1 

T = 8.D+0 

AA = 0.2 

Xy * -4.64 

y F = 2.2 

MMMMM = 23 
LOSC = 5 

SUM = 4.27375D-1 
SUMA = 7-5575D-1 
TIME = 9.9294D-2 
ANALYTICAL * 4.27375D-1 



177 



Case 86. 




SU8R 3 UTI ‘IS VALUE ( X , Y , T , A L , 8E ,tJ x V , , NC ALL , MC ALL ) 

THIS IS TiRLE ENTRY NUMBER 90 TP LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL-8 ( A,3,0 -h, 3-Z) 

IMPLICI T COMPLEX* lo ( C) 

DIMENSION NC ALL ( 5 ) * MCALL ( 5 I 

MMMMM = MMMMM + 1 

ZERO = 0.0+0 

CS1 = DCM=LX(X,Y) 

CS2 = OCMPLXIAL jZERO) 

C S3 = DCMPLXf 6E.ZER0) 

XPTR = X*T 
XPTI = Y»T 

CXPT = DCMPLXI XPTR, XPTI ) 

CS4 = CS1 + CS2 
CS5 = CS1 + C S3 
CS6 = CS4/CS5 

CALL CLOG (CS6, CST, I, mcaLL) 

CRESLT - CS7*C0SXP(CXPT) 

U = DREAL(CRESLT) 

V = OIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
BE = 2.5D-1 
T ■ 8.D+0 
AA * 0.2 
X p = -4.21 
y F = 2.1 
MMMMM =22 
LOSC = 5 

SUM = -2.90680D-2 
SUMA = 1.01743D-1 
TIME = 7.0096D-2 
ANALYTICAL = -2.90680D-2 
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Case 87 




In [(«* + 

2s 



Ci (at) 



SUB ? UT I^E VAL'JS (X, V, T f Ms U, V, MMMMM, MC ALL ) 

C THIS IS TABLE ENTRY NUM9E* 91 “>F LAPLACE TRANSFORMS BY SPIEGEL 

IMPLir I T PEAL-3 (AtS T D-H,n-Z) 

IM£>LICI t COMPLEX* 16 (C> 

0 1 ME NS I ON MC ALL (5) 

MMMMM = mmmmm + 1 
ZERO = 0 • D+O 
CS1 = OCMOLX(X,Y) 

CS2 - OCMP LX ( AL ,ZER3 ) 

CS3 = ((C*1*CS1) + (CS2*C$2))/( CS2*CS2) 

XP*P - X 
X°T I = Y *7 

CXPT = OC^PLX (XPT RtXPT I J 
CALL 3LOG ( CS3 » C$4, 1, *CALL) 

CNUM = CS4*CDEX?( CXPT ) 

CS5 = DC M PLX(2*0+0, ZERO) 

CRESLT = C NUM/ ( CS 5*C S 1 ) 

U = ORE AL ( CRESLT J 
V = DI M AG( CPESLTi 
RETJPN 
END 



AL * 1.25D-1 
T * 8. D+O 
AA = 0.2 
X F = -5.09 

y F * 2.3 

MMMMM = 2k 
LOSC = 5 

SUM * -3.37^0UD-1 
SUMA = 1.27603D+0 
TIME * 8.82UUD-2 
ANALYTICAL = -3O7^0^D-l 
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Case 88. 



I 




In {(* + a)/a\ 

9 



Ei (at) 



SUBROUTINE VALUE (X, Y, T, Alt U, V, MMMMM f MCALL) 
THIS IS TABLE ENTRY NUMBER 92 OF LAPLACE TRANSFORMS 
IMPLICIT REAL-8 ( A t 3 , 0-M, 0- Z ) 

I PPLI C I T COMPLEX* 16 ( C) 

DIMENSION ^CALL ( 5) 

MMMMM = MMMMM ♦ j. 

ZERO = C.D+O 
CS1 = DC M0 L X( Xf Y ) 

CS2 = DCMOLXt ALt ZERO) 

C S3 = ( CSl + CS2) /CS2 
XPTR = X*T 
X FTI = Y* T 

CXPT = DCMPLXtXPTRtXPTI J 
CALL CLOG (CS3t C S A » 1, MCALL) 

CNUM = C$A*CDtXP( CX°T ) 

CRESL T = CNUM /C SI 
U * DPEAL (C p ESLT) 

V = DIMAG(CRESLT) 

RETURN 

END 



AL = 1.25D-1 
T 8.D+0 
AA * 0.2 
x F = -5.09 
y F E 2.3 

MMMMM = 24 
LOSC = 5 

SUM - 2 . 193 B^D -1 
SUMA * 1.5509^0+0 
TIME * 7.212UD-2 
ANALYTICAL = 2.19384D-1 



BY SPIEGFL 
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Case 89 




- (In t + y) 

y = Euler’s constant = .5772156... 



SUBROUTINE value (X, Y, T, AL » U, V, mmmmm, MCALL) 

C THIS IS TABLE ENTRY NUMBER 96 CF LAPLACE TRANSFORMS BY SPIEGEL 

IMPLICIT 3 EM. *8 { A ,3 ,0-rl, 3-Z ) 

I MPL I C IT C 0 M P L EX * 1 6 (C) 

DIMENSION “CALL ( 5) 

MMMMM 3 MV VI MM + 1 

ZERO = O.O+O 
CS1 * OCMPLXI X f Y) 

CALL CLOG (CSU CS2, i» MCALL) 

X PTR = X*t 
XPTI = Y*T 

CXPT = OCMPLXIXPTPt XPTI ) 

CNUM = C. c 2*CDtXP( CXPT ) 

CRESLT = CNUM/CSI 
U = OR E AL ( C° ESLT ) 

V = DI MAG(CRESLT) 

RETURN 

ENO 



AL = 1.25D-1 
T = 8 . D+0 
AA = 0.2 
X F = -5.09 
y F = 2.3 

MMMMM = 24 
LOSC = 5 

SUM = -2.6566^D+0 
SUMA = 2.80606D+0 
TIME = 7.7818D-2 
ANALYTICAL = -2.656660+0 
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Case 90. 




(In I + 7 )= - Jr’ 

7 = Eulor’a conjtunl = .5772156... 



SUBROUTINE value (X, Y, T, SL, U, V, MMMmm, MCALU 



IMPLICI T REAL*8 ( A ,8 ,0-H , 0 - 1 J 
IMPL IC IT Cn M °LEX«=I6 <C ) 
DIMENSION MC ALL 151 

MMMMM s MMMMM + X 

ZERO = O.D+O 
CS1 = OCMPLXIX.Y) 

XPTR * X*T 
X PT I = Y*T 

CXPT = DC MP LX (X P T R , X P TI ) 

CALL CL P G (CS1, CS2, 1, M C A L L ) 
CMJM = CS2*CS2’ ! ‘C0EXP(CXP T ) 
CRESLT = CNUM/CSl 
U * DR EAL< CPESLT) 

V = D I M AG (CRESLT ) 

RETURN 

END 



C 



THIS IS TABLE ENTRY NUMBER '57 OF LAPLACE TRANSFORMS BY SPIEGEL 



AL = 1.25D-1 



T * 8.D+0 



AA “ 0.2 




MMMMM = 24 



LOSC = 4 



SUM = 5.U1289D+0 
SUMA = 5.43^?D+0 
TIME = 7.7376D-2 



ANALYTICAL = 5.41289D+0 
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Case 91 . 




SUBPOITINE VALUE (XtY,T f AL,BE,U,V t MMMMM,^CALLfMCALL) 
C THIS IS TABLE ENTRY *21 RAGE 265 ROBERTS ANO KAUFMAN 

IMPLICI T -EAL*8 { A,6,D-H, 1-L) 

I VPLICI " COMPLEX* Id ( C) 

DIMENSION MCALL(5)t MCALL15) 

MMMMM = + l 

ZERO = O.D + O 
CS1 * DCMC>LX( X,Y) 

CS2 = DCM& LX { AL t ZERO ) 

XP7R = X* T 
XPTI = Y* T 

CXPT = DCMPLX (XPTR f XPT I) 

C S3 = ( <CS1*CS1) ♦ ( C S2*C S2 ) ) 

CALL CPDWS* <CS3t C05N3M f 5.0-lt 1,NCALL) 

C$4 = (COENOm + CSD/CS2 
3ALL CLOG ( C S A t CS5t If *CALL) 

C NUM - CS5 *CDEX D ( CX^T ) 

CRESL T = CNUM/COENOM 
U = ORE AL { ESLT) 

V = 0IMAG{ CRESLT) 

RE T URN 

ENO 



AL = 1.25D-1 
T * 8.D+0 
AA = 0.3 
X F = -4.99 
y F - 2.3 
MMMMM = 24 
LOSC = 5 

SUM = -1.38634D-1 
SUMA = 3.92274D+0 
TIME = 1.2126D+0 
ANALYTICAL = -1.38634D-1 
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Case 92 



This 

page 





[- 




4 



[s* f « X J ^ 




f H„ M 



is known as the Struve function. See ref. (1), 
496. 



SUBROUTINE VALUE ( X , Y , T ,A L ,3 z , U , V . H“mmm , NC ALL , MC ALL ) 
C T HIS IS T ABLE 5N - RY *12 R*GE 264 ROBERTS AND KAUFMAN 

IMPLICIT R EAL*ii ( A ,3 , D-H , D-2) 

IMPLIC It C 0 MO L EX * 16 (C) 

DIMENSION MCALL(S), MCALL ( 5 ) 

MMMMM = MMMMM *. I 

ZERO = 0.0+0 
CS1 = OC MPLX (X » Y ) 

CS2 = DCMP LXt AL »ZER0) 

XPTR = X 
XPTI = Y** 

CXPT * OCMPLX(XRTR,XPTI ) 

CS3 = ((CS1*CSI) + ( C S 2 *CS2 ) ) 

CALL CPOWEF. (03 3 » COEND'* » 5. D-l , L * NCALL ) 

CS4 = < CDENOm + S S2 J /C SI 
CALL CL r G (CS4,CS5,i .MCALLI 

CNUM = CS5*G0cXP(CX?T) 

CRESLT = CNUM/CDENDM 
U « DREAL(C 3 ESL~J 
V = DIMAG(CRESLT) 

RETJRN 

END 



AL «= 1.25D-1 
T = 8.D+0 
AA * 0.3 
Xy = -k.5k 
y F = 2.2 
MMMMM = 23 
LOSC = 4 

SUM = 8.932440-1 
SUMA = 9.27960D-1 
TIME = 1.1341D-1 
ANALYTICAL = 8.93244D-1 
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Case 93 



Heaviside’s unit function — o) 




Fit) 



0 



SU8RDUT INE VALUE < X , Y , T , A L , U# V , mm MMM) 

THIS IS TABLE ENTRY NUMBER 139 PF LAPLACE TRANSFORMS BY SPIEGEL 
IMPLICIT REAL-8 (A,a.0-H,0-2) 

IMPLICIT CPMPLEX*16 ( C ) 

MMMMM = MMMMM * 1 
CS s DC MPL X ( X f Y ) 

X P T R = X * ( T — AL ) 

XPTI = Y*( T -AL) 

CXPT = DCMPLXtXPTRtXPTI ) 

CRESLT = CDEXPICXPT )/CS 
U = DRE AL(CRE^LT) 

V = DIMAG(CRESLT) 

retjrn 

END 



AL = 1.25D-1 
T = 8.0625D+0 
AA = 0.2 
X p = ~k.6k 

P 

y F = 2 * 2 
MMMMM = 23 

LOSC = 5 

SUM = 1 . 00000D+0 
SUMA = 1. ** 16820+0 
TIME = 6.055^0-2 
ANALYTICAL = l.OOOOOD+O 
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Case 94 




SUBROUTINE VALUE (X ,Y ,T , AL. U , V , YMMMm j 
C THIS TRANSFORM WAS TAKEN FROM R0BER t S ANO KAUFMAN 

IMPLICIT °EAL*8 (4-H.O-Z) 

MMMMM = MMMMM ♦ 1 

TEM 1 = -1.0+0«=AL«AL=»X 
T EM2 = AL*AL*Y 

TEM3 = A.D*C*< ( X«X) +• <Y*Y )) 

TEMA x TEM 1 /TEM 3 
TEM5 = T E M 2 / T EM3 
TEM6 = OEXPITEMAJxOCaSITSMF) 

TEM? = 0EX°(TEM4) *OSIN(TEM5) 

TEM8 = (TgM6*XI + C" E u 7 xy ) 

TEM9 * ( TE M 7* X) - ITSM6*Y) 

TEM10 = <X*X) * { Y*Y ) 

TE Ml 1 = TEM8/TEM10 

TEM 12 = TE M< 3/TE M 1 0 

TEM13 x OEXP(X*T)»OCOS(Y»’" J 

TEM1A x OEXP( X*T) =*D $ I N ( Y* T J 

U= ( TEM 13*TEM 1 1 J - ( TEMlAxTE M12) 

V = ( TE Ml A * T EMI 1 ) + ( T £M13» T E M 12 ) 

RETURN 

END 



AL = 2.5D-1 
T *= 1.6D+1 
AA = 0.2 
Xp = -3.04 
y F = 1.8 
MMMMM = 19 
LOSC = 8 

SUM = 7.65198D-1 
SUMA = 3.98573D+0 
TIME = 6 .36480-2 
ANALYTICAL = 7 . 65198 D-I 
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Case 95 





- at 

i - iatc 



SUBROUTINE VALUE (X ,Y ,T ,AL,U , V t 

THIS TRANSFORM WAS TAKEN FRO* ROBERTS AND KAUFMAN 
IMPLICIT REAL-3 I A, 3 , 0-H, 0-Z ) 

IMPLICI T COMPLEX* 16 ( C) 

MMMMM = mmmmm ♦ 1 
ZERO = 0 * D> 0 
CS1 = OC MP L X( X, Y) 

CS2 = DCMPLX ( AL, ZERO) 

C S3 = CS1 - CS2 
CS4 = CS1 + CS2 
XPTR * X* T 
XPTI = Y*T 

CXPT = DC.MPLX(XPTR t XPTI) 

CNUM = CS3*CS3*CJEXP( CXPT) 

COENOM = CS I*C S 4* CS 4 
C RES L~ = CNUM/CJENGM 
U = 0 RE AL ( C RESL T ) 

V = DI M AGt CRESLT) 

RET J RN 
END 



AL * 1.25D-1 
T = 8.D+0 
AA = 0.2 
Xp * -4.6U 

y F - 2.2 

MMMMM = 23 

Lose = 5 

SUM * -4.71518D-1 
SUMA = 5-0686^0-1 
TIME = 6.91B6D-2 
ANALYTICAL = -4.715183-1 
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